!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!

! **************************************************************************************************
!> \par History
!>      - 02.2004 flexible normalization of basis sets [jgh]
!>      - 07.2014 Add a set of contraction coefficient that only work on active
!>                functions
!> \author Matthias Krack (04.07.2000)
! **************************************************************************************************
MODULE basis_set_types

   USE ai_coulomb,                      ONLY: coulomb2
   USE bibliography,                    ONLY: VandeVondele2007,&
                                              cite_reference
   USE cp_linked_list_input,            ONLY: cp_sll_val_next,&
                                              cp_sll_val_type
   USE cp_parser_methods,               ONLY: parser_get_object,&
                                              parser_search_string
   USE cp_parser_types,                 ONLY: cp_parser_type,&
                                              parser_create,&
                                              parser_release
   USE input_section_types,             ONLY: section_vals_list_get,&
                                              section_vals_type,&
                                              section_vals_val_get
   USE input_val_types,                 ONLY: val_get,&
                                              val_type
   USE kinds,                           ONLY: default_path_length,&
                                              default_string_length,&
                                              dp
   USE mathconstants,                   ONLY: dfac,&
                                              pi
   USE memory_utilities,                ONLY: reallocate
   USE message_passing,                 ONLY: mp_para_env_type
   USE orbital_pointers,                ONLY: coset,&
                                              indco,&
                                              init_orbital_pointers,&
                                              nco,&
                                              ncoset,&
                                              nso,&
                                              nsoset
   USE orbital_symbols,                 ONLY: cgf_symbol,&
                                              sgf_symbol
   USE orbital_transformation_matrices, ONLY: init_spherical_harmonics,&
                                              orbtramat
   USE sto_ng,                          ONLY: get_sto_ng
   USE string_utilities,                ONLY: integer_to_string,&
                                              remove_word,&
                                              uppercase
   USE util,                            ONLY: sort
#include "../base/base_uses.f90"

   IMPLICIT NONE

   PRIVATE

   ! Global parameters (only in this module)

   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'basis_set_types'

   ! basis set sort criteria
   INTEGER, PARAMETER, PUBLIC               :: basis_sort_default = 0, &
                                               basis_sort_zet = 1

! **************************************************************************************************
   ! Define the Gaussian-type orbital basis set type

   TYPE gto_basis_set_type
      !MK PRIVATE
      CHARACTER(LEN=default_string_length)       :: name = ""
      CHARACTER(LEN=default_string_length)       :: aliases = ""
      REAL(KIND=dp)                              :: kind_radius = 0.0_dp
      REAL(KIND=dp)                              :: short_kind_radius = 0.0_dp
      INTEGER                                    :: norm_type = -1
      INTEGER                                    :: ncgf = -1, nset = -1, nsgf = -1
      CHARACTER(LEN=12), DIMENSION(:), POINTER   :: cgf_symbol => NULL()
      CHARACTER(LEN=6), DIMENSION(:), POINTER    :: sgf_symbol => NULL()
      REAL(KIND=dp), DIMENSION(:), POINTER       :: norm_cgf => NULL(), set_radius => NULL()
      INTEGER, DIMENSION(:), POINTER             :: lmax => NULL(), lmin => NULL(), &
                                                    lx => NULL(), ly => NULL(), lz => NULL(), &
                                                    m => NULL(), ncgf_set => NULL(), &
                                                    npgf => NULL(), nsgf_set => NULL(), nshell => NULL()
      REAL(KIND=dp), DIMENSION(:, :), POINTER    :: cphi => NULL(), pgf_radius => NULL(), sphi => NULL(), &
                                                    scon => NULL(), zet => NULL()
      INTEGER, DIMENSION(:, :), POINTER          :: first_cgf => NULL(), first_sgf => NULL(), l => NULL(), &
                                                    last_cgf => NULL(), last_sgf => NULL(), n => NULL()
      REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc => NULL()
   END TYPE gto_basis_set_type

   TYPE gto_basis_set_p_type
      TYPE(gto_basis_set_type), POINTER :: gto_basis_set => NULL()
   END TYPE gto_basis_set_p_type

! **************************************************************************************************
   ! Define the Slater-type orbital basis set type

   TYPE sto_basis_set_type
      PRIVATE
      CHARACTER(LEN=default_string_length)       :: name = ""
      INTEGER                                    :: nshell = -1
      CHARACTER(LEN=6), DIMENSION(:), POINTER    :: symbol => NULL()
      INTEGER, DIMENSION(:), POINTER             :: nq => NULL(), lq => NULL()
      REAL(KIND=dp), DIMENSION(:), POINTER       :: zet => NULL()
   END TYPE sto_basis_set_type

! **************************************************************************************************
   INTERFACE read_gto_basis_set
      MODULE PROCEDURE read_gto_basis_set1, read_gto_basis_set2
   END INTERFACE
! **************************************************************************************************

   ! Public subroutines
   PUBLIC :: allocate_gto_basis_set, &
             deallocate_gto_basis_set, &
             get_gto_basis_set, &
             init_aux_basis_set, &
             init_cphi_and_sphi, &
             init_orb_basis_set, &
             read_gto_basis_set, &
             copy_gto_basis_set, &
             create_primitive_basis_set, &
             combine_basis_sets, &
             set_gto_basis_set, &
             sort_gto_basis_set, &
             write_gto_basis_set, &
             dump_gto_basis_set, &
             write_orb_basis_set

   PUBLIC :: allocate_sto_basis_set, &
             read_sto_basis_set, &
             create_gto_from_sto_basis, &
             deallocate_sto_basis_set, &
             set_sto_basis_set, &
             srules, process_gto_basis

   ! Public data types
   PUBLIC :: gto_basis_set_p_type, &
             gto_basis_set_type, &
             sto_basis_set_type

CONTAINS

! **************************************************************************************************
!> \brief ...
!> \param gto_basis_set ...
! **************************************************************************************************
   SUBROUTINE allocate_gto_basis_set(gto_basis_set)

      ! Allocate a Gaussian-type orbital (GTO) basis set data set.

      ! - Creation (26.10.2000,MK)

      TYPE(gto_basis_set_type), POINTER                  :: gto_basis_set

      CALL deallocate_gto_basis_set(gto_basis_set)

      ALLOCATE (gto_basis_set)

   END SUBROUTINE allocate_gto_basis_set

! **************************************************************************************************
!> \brief ...
!> \param gto_basis_set ...
! **************************************************************************************************
   SUBROUTINE deallocate_gto_basis_set(gto_basis_set)

      ! Deallocate a Gaussian-type orbital (GTO) basis set data set.

      ! - Creation (03.11.2000,MK)

      TYPE(gto_basis_set_type), POINTER                  :: gto_basis_set

      IF (ASSOCIATED(gto_basis_set)) THEN
         IF (ASSOCIATED(gto_basis_set%cgf_symbol)) DEALLOCATE (gto_basis_set%cgf_symbol)
         IF (ASSOCIATED(gto_basis_set%sgf_symbol)) DEALLOCATE (gto_basis_set%sgf_symbol)
         IF (ASSOCIATED(gto_basis_set%norm_cgf)) DEALLOCATE (gto_basis_set%norm_cgf)
         IF (ASSOCIATED(gto_basis_set%set_radius)) DEALLOCATE (gto_basis_set%set_radius)
         IF (ASSOCIATED(gto_basis_set%lmax)) DEALLOCATE (gto_basis_set%lmax)
         IF (ASSOCIATED(gto_basis_set%lmin)) DEALLOCATE (gto_basis_set%lmin)
         IF (ASSOCIATED(gto_basis_set%lx)) DEALLOCATE (gto_basis_set%lx)
         IF (ASSOCIATED(gto_basis_set%ly)) DEALLOCATE (gto_basis_set%ly)
         IF (ASSOCIATED(gto_basis_set%lz)) DEALLOCATE (gto_basis_set%lz)
         IF (ASSOCIATED(gto_basis_set%m)) DEALLOCATE (gto_basis_set%m)
         IF (ASSOCIATED(gto_basis_set%ncgf_set)) DEALLOCATE (gto_basis_set%ncgf_set)
         IF (ASSOCIATED(gto_basis_set%npgf)) DEALLOCATE (gto_basis_set%npgf)
         IF (ASSOCIATED(gto_basis_set%nsgf_set)) DEALLOCATE (gto_basis_set%nsgf_set)
         IF (ASSOCIATED(gto_basis_set%nshell)) DEALLOCATE (gto_basis_set%nshell)
         IF (ASSOCIATED(gto_basis_set%cphi)) DEALLOCATE (gto_basis_set%cphi)
         IF (ASSOCIATED(gto_basis_set%pgf_radius)) DEALLOCATE (gto_basis_set%pgf_radius)
         IF (ASSOCIATED(gto_basis_set%sphi)) DEALLOCATE (gto_basis_set%sphi)
         IF (ASSOCIATED(gto_basis_set%scon)) DEALLOCATE (gto_basis_set%scon)
         IF (ASSOCIATED(gto_basis_set%zet)) DEALLOCATE (gto_basis_set%zet)
         IF (ASSOCIATED(gto_basis_set%first_cgf)) DEALLOCATE (gto_basis_set%first_cgf)
         IF (ASSOCIATED(gto_basis_set%first_sgf)) DEALLOCATE (gto_basis_set%first_sgf)
         IF (ASSOCIATED(gto_basis_set%l)) DEALLOCATE (gto_basis_set%l)
         IF (ASSOCIATED(gto_basis_set%last_cgf)) DEALLOCATE (gto_basis_set%last_cgf)
         IF (ASSOCIATED(gto_basis_set%last_sgf)) DEALLOCATE (gto_basis_set%last_sgf)
         IF (ASSOCIATED(gto_basis_set%n)) DEALLOCATE (gto_basis_set%n)
         IF (ASSOCIATED(gto_basis_set%gcc)) DEALLOCATE (gto_basis_set%gcc)
         DEALLOCATE (gto_basis_set)
      END IF
   END SUBROUTINE deallocate_gto_basis_set

! **************************************************************************************************
!> \brief ...
!> \param basis_set_in ...
!> \param basis_set_out ...
! **************************************************************************************************
   SUBROUTINE copy_gto_basis_set(basis_set_in, basis_set_out)

      ! Copy a Gaussian-type orbital (GTO) basis set data set.

      TYPE(gto_basis_set_type), INTENT(IN)               :: basis_set_in
      TYPE(gto_basis_set_type), POINTER                  :: basis_set_out

      INTEGER                                            :: maxco, maxpgf, maxshell, ncgf, nset, nsgf

      CALL allocate_gto_basis_set(basis_set_out)

      basis_set_out%name = basis_set_in%name
      basis_set_out%aliases = basis_set_in%aliases
      basis_set_out%kind_radius = basis_set_in%kind_radius
      basis_set_out%norm_type = basis_set_in%norm_type
      basis_set_out%nset = basis_set_in%nset
      basis_set_out%ncgf = basis_set_in%ncgf
      basis_set_out%nsgf = basis_set_in%nsgf
      nset = basis_set_in%nset
      ncgf = basis_set_in%ncgf
      nsgf = basis_set_in%nsgf
      ALLOCATE (basis_set_out%cgf_symbol(ncgf))
      ALLOCATE (basis_set_out%sgf_symbol(nsgf))
      basis_set_out%cgf_symbol = basis_set_in%cgf_symbol
      basis_set_out%sgf_symbol = basis_set_in%sgf_symbol
      ALLOCATE (basis_set_out%norm_cgf(ncgf))
      basis_set_out%norm_cgf = basis_set_in%norm_cgf
      ALLOCATE (basis_set_out%set_radius(nset))
      basis_set_out%set_radius = basis_set_in%set_radius
      ALLOCATE (basis_set_out%lmax(nset), basis_set_out%lmin(nset), basis_set_out%npgf(nset), basis_set_out%nshell(nset))
      basis_set_out%lmax = basis_set_in%lmax
      basis_set_out%lmin = basis_set_in%lmin
      basis_set_out%npgf = basis_set_in%npgf
      basis_set_out%nshell = basis_set_in%nshell
      ALLOCATE (basis_set_out%lx(ncgf), basis_set_out%ly(ncgf), basis_set_out%lz(ncgf), basis_set_out%m(nsgf))
      basis_set_out%lx = basis_set_in%lx
      basis_set_out%ly = basis_set_in%ly
      basis_set_out%lz = basis_set_in%lz
      basis_set_out%m = basis_set_in%m
      ALLOCATE (basis_set_out%ncgf_set(nset), basis_set_out%nsgf_set(nset))
      basis_set_out%ncgf_set = basis_set_in%ncgf_set
      basis_set_out%nsgf_set = basis_set_in%nsgf_set
      maxco = SIZE(basis_set_in%cphi, 1)
      ALLOCATE (basis_set_out%cphi(maxco, ncgf), basis_set_out%sphi(maxco, nsgf), basis_set_out%scon(maxco, nsgf))
      basis_set_out%cphi = basis_set_in%cphi
      basis_set_out%sphi = basis_set_in%sphi
      basis_set_out%scon = basis_set_in%scon
      maxpgf = MAXVAL(basis_set_in%npgf)
      ALLOCATE (basis_set_out%pgf_radius(maxpgf, nset), basis_set_out%zet(maxpgf, nset))
      basis_set_out%pgf_radius = basis_set_in%pgf_radius
      basis_set_out%zet = basis_set_in%zet
      maxshell = MAXVAL(basis_set_in%nshell)
      ALLOCATE (basis_set_out%first_cgf(maxshell, nset), basis_set_out%first_sgf(maxshell, nset))
      ALLOCATE (basis_set_out%last_cgf(maxshell, nset), basis_set_out%last_sgf(maxshell, nset))
      basis_set_out%first_cgf = basis_set_in%first_cgf
      basis_set_out%first_sgf = basis_set_in%first_sgf
      basis_set_out%last_cgf = basis_set_in%last_cgf
      basis_set_out%last_sgf = basis_set_in%last_sgf
      ALLOCATE (basis_set_out%n(maxshell, nset), basis_set_out%l(maxshell, nset))
      basis_set_out%n = basis_set_in%n
      basis_set_out%l = basis_set_in%l
      ALLOCATE (basis_set_out%gcc(maxpgf, maxshell, nset))
      basis_set_out%gcc = basis_set_in%gcc

   END SUBROUTINE copy_gto_basis_set

! **************************************************************************************************
!> \brief ...
!> \param basis_set ...
!> \param pbasis ...
!> \param lmax ...
! **************************************************************************************************
   SUBROUTINE create_primitive_basis_set(basis_set, pbasis, lmax)

      ! Create a primitives only basis set

      TYPE(gto_basis_set_type), INTENT(IN)               :: basis_set
      TYPE(gto_basis_set_type), POINTER                  :: pbasis
      INTEGER, INTENT(IN), OPTIONAL                      :: lmax

      INTEGER                                            :: i, ico, ip, ipgf, iset, ishell, l, lm, &
                                                            lshell, m, maxco, mpgf, nc, ncgf, ns, &
                                                            nset, nsgf
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: nindex, nprim
      REAL(KIND=dp)                                      :: zet0
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: zet, zeta

      mpgf = SUM(basis_set%npgf)
      lm = MAXVAL(basis_set%lmax)
      ALLOCATE (zet(mpgf, 0:lm), zeta(mpgf, lm + 1), nindex(mpgf), nprim(0:lm))
      zet = 0.0_dp
      zeta = 0.0_dp
      DO l = 0, lm
         ip = 0
         DO iset = 1, basis_set%nset
            IF (basis_set%lmin(iset) <= l .AND. basis_set%lmax(iset) >= l) THEN
               DO ipgf = 1, basis_set%npgf(iset)
                  ip = ip + 1
                  zet(ip, l) = basis_set%zet(ipgf, iset)
               END DO
            END IF
         END DO
         nprim(l) = ip
      END DO

      ! sort exponents
      DO l = 0, lm
         zet(1:nprim(l), l) = -zet(1:nprim(l), l)
         CALL sort(zet(1:nprim(l), l), nprim(l), nindex)
         ! remove duplicates
         ip = 0
         zet0 = 0.0_dp
         DO i = 1, nprim(l)
            IF (ABS(zet0 - zet(i, l)) > 1.e-6_dp) THEN
               ip = ip + 1
               zeta(ip, l + 1) = zet(i, l)
            END IF
         END DO
         nprim(l) = ip
         !
         zeta(1:ip, l + 1) = -zeta(1:ip, l + 1)
      END DO

      CALL allocate_gto_basis_set(pbasis)

      IF (PRESENT(lmax)) THEN
         ! if requested, reduce max l val
         lm = MIN(lm, lmax)
      END IF

      IF (LEN_TRIM(basis_set%name) + 10 > default_string_length) THEN
         CPWARN("The name of the primitive basis set will be truncated.")
      END IF
      pbasis%name = TRIM(basis_set%name)//"_primitive"
      pbasis%kind_radius = basis_set%kind_radius
      pbasis%short_kind_radius = basis_set%short_kind_radius
      pbasis%norm_type = basis_set%norm_type
      nset = lm + 1
      pbasis%nset = nset
      ALLOCATE (pbasis%lmax(nset), pbasis%lmin(nset), pbasis%npgf(nset), pbasis%nshell(nset))
      DO iset = 1, nset
         pbasis%lmax(iset) = iset - 1
         pbasis%lmin(iset) = iset - 1
         pbasis%npgf(iset) = nprim(iset - 1)
         pbasis%nshell(iset) = nprim(iset - 1)
      END DO
      pbasis%ncgf = 0
      pbasis%nsgf = 0
      DO l = 0, lm
         pbasis%ncgf = pbasis%ncgf + nprim(l)*((l + 1)*(l + 2))/2
         pbasis%nsgf = pbasis%nsgf + nprim(l)*(2*l + 1)
      END DO
      mpgf = MAXVAL(nprim)
      ALLOCATE (pbasis%zet(mpgf, nset))
      pbasis%zet(1:mpgf, 1:nset) = zeta(1:mpgf, 1:nset)

      ALLOCATE (pbasis%l(mpgf, nset), pbasis%n(mpgf, nset))
      DO iset = 1, nset
         DO ip = 1, nprim(iset - 1)
            pbasis%l(ip, iset) = iset - 1
            pbasis%n(ip, iset) = iset + ip - 1
         END DO
      END DO

      ALLOCATE (pbasis%cgf_symbol(pbasis%ncgf))
      ALLOCATE (pbasis%lx(pbasis%ncgf))
      ALLOCATE (pbasis%ly(pbasis%ncgf))
      ALLOCATE (pbasis%lz(pbasis%ncgf))
      ALLOCATE (pbasis%m(pbasis%nsgf))
      ALLOCATE (pbasis%sgf_symbol(pbasis%nsgf))
      ALLOCATE (pbasis%ncgf_set(nset), pbasis%nsgf_set(nset))

      ncgf = 0
      nsgf = 0
      DO iset = 1, nset
         l = iset - 1
         pbasis%ncgf_set(iset) = nprim(l)*((l + 1)*(l + 2))/2
         pbasis%nsgf_set(iset) = nprim(l)*(2*l + 1)
         DO ishell = 1, pbasis%nshell(iset)
            lshell = pbasis%l(ishell, iset)
            DO ico = ncoset(lshell - 1) + 1, ncoset(lshell)
               ncgf = ncgf + 1
               pbasis%lx(ncgf) = indco(1, ico)
               pbasis%ly(ncgf) = indco(2, ico)
               pbasis%lz(ncgf) = indco(3, ico)
               pbasis%cgf_symbol(ncgf) = &
                  cgf_symbol(pbasis%n(ishell, iset), [pbasis%lx(ncgf), pbasis%ly(ncgf), pbasis%lz(ncgf)])
            END DO
            DO m = -lshell, lshell
               nsgf = nsgf + 1
               pbasis%m(nsgf) = m
               pbasis%sgf_symbol(nsgf) = sgf_symbol(pbasis%n(ishell, iset), lshell, m)
            END DO
         END DO
      END DO
      CPASSERT(ncgf == pbasis%ncgf)
      CPASSERT(nsgf == pbasis%nsgf)

      ALLOCATE (pbasis%gcc(mpgf, mpgf, nset))
      pbasis%gcc = 0.0_dp
      DO iset = 1, nset
         DO i = 1, mpgf
            pbasis%gcc(i, i, iset) = 1.0_dp
         END DO
      END DO

      ALLOCATE (pbasis%first_cgf(mpgf, nset))
      ALLOCATE (pbasis%first_sgf(mpgf, nset))
      ALLOCATE (pbasis%last_cgf(mpgf, nset))
      ALLOCATE (pbasis%last_sgf(mpgf, nset))
      nc = 0
      ns = 0
      maxco = 0
      DO iset = 1, nset
         DO ishell = 1, pbasis%nshell(iset)
            lshell = pbasis%l(ishell, iset)
            pbasis%first_cgf(ishell, iset) = nc + 1
            nc = nc + nco(lshell)
            pbasis%last_cgf(ishell, iset) = nc
            pbasis%first_sgf(ishell, iset) = ns + 1
            ns = ns + nso(lshell)
            pbasis%last_sgf(ishell, iset) = ns
         END DO
         maxco = MAX(maxco, pbasis%npgf(iset)*ncoset(pbasis%lmax(iset)))
      END DO

      ALLOCATE (pbasis%norm_cgf(ncgf))
      ALLOCATE (pbasis%cphi(maxco, ncgf))
      pbasis%cphi = 0.0_dp
      ALLOCATE (pbasis%sphi(maxco, nsgf))
      pbasis%sphi = 0.0_dp
      ALLOCATE (pbasis%scon(maxco, ncgf))
      pbasis%scon = 0.0_dp
      ALLOCATE (pbasis%set_radius(nset))
      ALLOCATE (pbasis%pgf_radius(mpgf, nset))
      pbasis%pgf_radius = 0.0_dp

      CALL init_orb_basis_set(pbasis)

      DEALLOCATE (zet, zeta, nindex, nprim)

   END SUBROUTINE create_primitive_basis_set

! **************************************************************************************************
!> \brief ...
!> \param basis_set ...
!> \param basis_set_add ...
! **************************************************************************************************
   SUBROUTINE combine_basis_sets(basis_set, basis_set_add)

      ! Combine two Gaussian-type orbital (GTO) basis sets.

      TYPE(gto_basis_set_type), INTENT(INOUT)            :: basis_set
      TYPE(gto_basis_set_type), INTENT(IN)               :: basis_set_add

      CHARACTER(LEN=12), DIMENSION(:), POINTER           :: cgf_symbol
      CHARACTER(LEN=6), DIMENSION(:), POINTER            :: sgf_symbol
      INTEGER                                            :: iset, ishell, lshell, maxco, maxpgf, &
                                                            maxshell, nc, ncgf, ncgfn, ncgfo, ns, &
                                                            nset, nsetn, nseto, nsgf, nsgfn, nsgfo

      IF (LEN_TRIM(basis_set%name) + LEN_TRIM(basis_set_add%name) > default_string_length) THEN
         CPWARN("The name of the combined GTO basis set will be truncated.")
      END IF
      basis_set%name = TRIM(basis_set%name)//TRIM(basis_set_add%name)
      basis_set%nset = basis_set%nset + basis_set_add%nset
      basis_set%ncgf = basis_set%ncgf + basis_set_add%ncgf
      basis_set%nsgf = basis_set%nsgf + basis_set_add%nsgf
      nset = basis_set%nset
      ncgf = basis_set%ncgf
      nsgf = basis_set%nsgf

      nsetn = basis_set_add%nset
      nseto = nset - nsetn
      CALL reallocate(basis_set%set_radius, 1, nset) ! to be defined later
      CALL reallocate(basis_set%lmax, 1, nset)
      CALL reallocate(basis_set%lmin, 1, nset)
      CALL reallocate(basis_set%npgf, 1, nset)
      CALL reallocate(basis_set%nshell, 1, nset)
      basis_set%lmax(nseto + 1:nset) = basis_set_add%lmax(1:nsetn)
      basis_set%lmin(nseto + 1:nset) = basis_set_add%lmin(1:nsetn)
      basis_set%npgf(nseto + 1:nset) = basis_set_add%npgf(1:nsetn)
      basis_set%nshell(nseto + 1:nset) = basis_set_add%nshell(1:nsetn)
      CALL reallocate(basis_set%ncgf_set, 1, nset)
      CALL reallocate(basis_set%nsgf_set, 1, nset)
      basis_set%ncgf_set(nseto + 1:nset) = basis_set_add%ncgf_set(1:nsetn)
      basis_set%nsgf_set(nseto + 1:nset) = basis_set_add%nsgf_set(1:nsetn)

      nsgfn = basis_set_add%nsgf
      nsgfo = nsgf - nsgfn
      ncgfn = basis_set_add%ncgf
      ncgfo = ncgf - ncgfn

      ALLOCATE (cgf_symbol(ncgf), sgf_symbol(nsgf))
      cgf_symbol(1:ncgfo) = basis_set%cgf_symbol(1:ncgfo)
      cgf_symbol(ncgfo + 1:ncgf) = basis_set_add%cgf_symbol(1:ncgfn)
      sgf_symbol(1:nsgfo) = basis_set%sgf_symbol(1:nsgfo)
      sgf_symbol(nsgfo + 1:nsgf) = basis_set_add%sgf_symbol(1:nsgfn)
      DEALLOCATE (basis_set%cgf_symbol, basis_set%sgf_symbol)
      ALLOCATE (basis_set%cgf_symbol(ncgf), basis_set%sgf_symbol(nsgf))
      basis_set%cgf_symbol = cgf_symbol
      basis_set%sgf_symbol = sgf_symbol
      DEALLOCATE (cgf_symbol, sgf_symbol)

      CALL reallocate(basis_set%lx, 1, ncgf)
      CALL reallocate(basis_set%ly, 1, ncgf)
      CALL reallocate(basis_set%lz, 1, ncgf)
      CALL reallocate(basis_set%m, 1, nsgf)
      basis_set%lx(ncgfo + 1:ncgf) = basis_set_add%lx(1:ncgfn)
      basis_set%ly(ncgfo + 1:ncgf) = basis_set_add%ly(1:ncgfn)
      basis_set%lz(ncgfo + 1:ncgf) = basis_set_add%lz(1:ncgfn)
      basis_set%m(nsgfo + 1:nsgf) = basis_set_add%m(1:nsgfn)

      maxpgf = MAXVAL(basis_set%npgf)
      CALL reallocate(basis_set%zet, 1, maxpgf, 1, nset)
      nc = SIZE(basis_set_add%zet, 1)
      DO iset = 1, nsetn
         basis_set%zet(1:nc, nseto + iset) = basis_set_add%zet(1:nc, iset)
      END DO

      maxshell = MAXVAL(basis_set%nshell)
      CALL reallocate(basis_set%l, 1, maxshell, 1, nset)
      CALL reallocate(basis_set%n, 1, maxshell, 1, nset)
      nc = SIZE(basis_set_add%l, 1)
      DO iset = 1, nsetn
         basis_set%l(1:nc, nseto + iset) = basis_set_add%l(1:nc, iset)
         basis_set%n(1:nc, nseto + iset) = basis_set_add%n(1:nc, iset)
      END DO

      CALL reallocate(basis_set%first_cgf, 1, maxshell, 1, nset)
      CALL reallocate(basis_set%first_sgf, 1, maxshell, 1, nset)
      CALL reallocate(basis_set%last_cgf, 1, maxshell, 1, nset)
      CALL reallocate(basis_set%last_sgf, 1, maxshell, 1, nset)
      nc = 0
      ns = 0
      DO iset = 1, nset
         DO ishell = 1, basis_set%nshell(iset)
            lshell = basis_set%l(ishell, iset)
            basis_set%first_cgf(ishell, iset) = nc + 1
            nc = nc + nco(lshell)
            basis_set%last_cgf(ishell, iset) = nc
            basis_set%first_sgf(ishell, iset) = ns + 1
            ns = ns + nso(lshell)
            basis_set%last_sgf(ishell, iset) = ns
         END DO
      END DO

      CALL reallocate(basis_set%gcc, 1, maxpgf, 1, maxshell, 1, nset)
      nc = SIZE(basis_set_add%gcc, 1)
      ns = SIZE(basis_set_add%gcc, 2)
      DO iset = 1, nsetn
         basis_set%gcc(1:nc, 1:ns, nseto + iset) = basis_set_add%gcc(1:nc, 1:ns, iset)
      END DO

      ! these arrays are determined later using initialization calls
      CALL reallocate(basis_set%norm_cgf, 1, ncgf)
      maxco = MAX(SIZE(basis_set%cphi, 1), SIZE(basis_set_add%cphi, 1))
      CALL reallocate(basis_set%cphi, 1, maxco, 1, ncgf)
      CALL reallocate(basis_set%sphi, 1, maxco, 1, nsgf)
      CALL reallocate(basis_set%scon, 1, maxco, 1, nsgf)
      CALL reallocate(basis_set%pgf_radius, 1, maxpgf, 1, nset)

   END SUBROUTINE combine_basis_sets

! **************************************************************************************************
!> \brief ...
!> \param gto_basis_set ...
!> \param name ...
!> \param aliases ...
!> \param norm_type ...
!> \param kind_radius ...
!> \param ncgf ...
!> \param nset ...
!> \param nsgf ...
!> \param cgf_symbol ...
!> \param sgf_symbol ...
!> \param norm_cgf ...
!> \param set_radius ...
!> \param lmax ...
!> \param lmin ...
!> \param lx ...
!> \param ly ...
!> \param lz ...
!> \param m ...
!> \param ncgf_set ...
!> \param npgf ...
!> \param nsgf_set ...
!> \param nshell ...
!> \param cphi ...
!> \param pgf_radius ...
!> \param sphi ...
!> \param scon ...
!> \param zet ...
!> \param first_cgf ...
!> \param first_sgf ...
!> \param l ...
!> \param last_cgf ...
!> \param last_sgf ...
!> \param n ...
!> \param gcc ...
!> \param maxco ...
!> \param maxl ...
!> \param maxpgf ...
!> \param maxsgf_set ...
!> \param maxshell ...
!> \param maxso ...
!> \param nco_sum ...
!> \param npgf_sum ...
!> \param nshell_sum ...
!> \param maxder ...
!> \param short_kind_radius ...
!> \param npgf_seg_sum number of primitives in "segmented contraction format"
! **************************************************************************************************
   SUBROUTINE get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, &
                                nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, &
                                m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, &
                                last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, &
                                npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum)

      ! Get informations about a Gaussian-type orbital (GTO) basis set.

      ! - Creation (10.01.2002,MK)

      TYPE(gto_basis_set_type), INTENT(IN)               :: gto_basis_set
      CHARACTER(LEN=default_string_length), &
         INTENT(OUT), OPTIONAL                           :: name, aliases
      INTEGER, INTENT(OUT), OPTIONAL                     :: norm_type
      REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: kind_radius
      INTEGER, INTENT(OUT), OPTIONAL                     :: ncgf, nset, nsgf
      CHARACTER(LEN=12), DIMENSION(:), OPTIONAL, POINTER :: cgf_symbol
      CHARACTER(LEN=6), DIMENSION(:), OPTIONAL, POINTER  :: sgf_symbol
      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: norm_cgf, set_radius
      INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: lmax, lmin, lx, ly, lz, m, ncgf_set, &
                                                            npgf, nsgf_set, nshell
      REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: cphi, pgf_radius, sphi, scon, zet
      INTEGER, DIMENSION(:, :), OPTIONAL, POINTER        :: first_cgf, first_sgf, l, last_cgf, &
                                                            last_sgf, n
      REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
         POINTER                                         :: gcc
      INTEGER, INTENT(OUT), OPTIONAL                     :: maxco, maxl, maxpgf, maxsgf_set, &
                                                            maxshell, maxso, nco_sum, npgf_sum, &
                                                            nshell_sum
      INTEGER, INTENT(IN), OPTIONAL                      :: maxder
      REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: short_kind_radius
      INTEGER, INTENT(OUT), OPTIONAL                     :: npgf_seg_sum

      INTEGER                                            :: iset, nder

      IF (PRESENT(name)) name = gto_basis_set%name
      IF (PRESENT(aliases)) aliases = gto_basis_set%aliases
      IF (PRESENT(norm_type)) norm_type = gto_basis_set%norm_type
      IF (PRESENT(kind_radius)) kind_radius = gto_basis_set%kind_radius
      IF (PRESENT(short_kind_radius)) short_kind_radius = gto_basis_set%short_kind_radius
      IF (PRESENT(ncgf)) ncgf = gto_basis_set%ncgf
      IF (PRESENT(nset)) nset = gto_basis_set%nset
      IF (PRESENT(nsgf)) nsgf = gto_basis_set%nsgf
      IF (PRESENT(cgf_symbol)) cgf_symbol => gto_basis_set%cgf_symbol
      IF (PRESENT(sgf_symbol)) sgf_symbol => gto_basis_set%sgf_symbol
      IF (PRESENT(norm_cgf)) norm_cgf => gto_basis_set%norm_cgf
      IF (PRESENT(set_radius)) set_radius => gto_basis_set%set_radius
      IF (PRESENT(lmax)) lmax => gto_basis_set%lmax
      IF (PRESENT(lmin)) lmin => gto_basis_set%lmin
      IF (PRESENT(lx)) lx => gto_basis_set%lx
      IF (PRESENT(ly)) ly => gto_basis_set%ly
      IF (PRESENT(lz)) lz => gto_basis_set%lz
      IF (PRESENT(m)) m => gto_basis_set%m
      IF (PRESENT(ncgf_set)) ncgf_set => gto_basis_set%ncgf_set
      IF (PRESENT(npgf)) npgf => gto_basis_set%npgf
      IF (PRESENT(nsgf_set)) nsgf_set => gto_basis_set%nsgf_set
      IF (PRESENT(nshell)) nshell => gto_basis_set%nshell
      IF (PRESENT(cphi)) cphi => gto_basis_set%cphi
      IF (PRESENT(pgf_radius)) pgf_radius => gto_basis_set%pgf_radius
      IF (PRESENT(sphi)) sphi => gto_basis_set%sphi
      IF (PRESENT(scon)) scon => gto_basis_set%scon
      IF (PRESENT(zet)) zet => gto_basis_set%zet
      IF (PRESENT(first_cgf)) first_cgf => gto_basis_set%first_cgf
      IF (PRESENT(first_sgf)) first_sgf => gto_basis_set%first_sgf
      IF (PRESENT(l)) l => gto_basis_set%l
      IF (PRESENT(last_cgf)) last_cgf => gto_basis_set%last_cgf
      IF (PRESENT(last_sgf)) last_sgf => gto_basis_set%last_sgf
      IF (PRESENT(n)) n => gto_basis_set%n
      IF (PRESENT(gcc)) gcc => gto_basis_set%gcc
      IF (PRESENT(maxco)) THEN
         maxco = 0
         IF (PRESENT(maxder)) THEN
            nder = maxder
         ELSE
            nder = 0
         END IF
         DO iset = 1, gto_basis_set%nset
            maxco = MAX(maxco, gto_basis_set%npgf(iset)* &
                        ncoset(gto_basis_set%lmax(iset) + nder))
         END DO
      END IF
      IF (PRESENT(maxl)) THEN
         maxl = -1
         DO iset = 1, gto_basis_set%nset
            maxl = MAX(maxl, gto_basis_set%lmax(iset))
         END DO
      END IF
      IF (PRESENT(maxpgf)) THEN
         maxpgf = 0
         DO iset = 1, gto_basis_set%nset
            maxpgf = MAX(maxpgf, gto_basis_set%npgf(iset))
         END DO
      END IF
      IF (PRESENT(maxsgf_set)) THEN
         maxsgf_set = 0
         DO iset = 1, gto_basis_set%nset
            maxsgf_set = MAX(maxsgf_set, gto_basis_set%nsgf_set(iset))
         END DO
      END IF
      IF (PRESENT(maxshell)) THEN ! MAXVAL on structure component avoided
         maxshell = 0
         DO iset = 1, gto_basis_set%nset
            maxshell = MAX(maxshell, gto_basis_set%nshell(iset))
         END DO
      END IF
      IF (PRESENT(maxso)) THEN
         maxso = 0
         DO iset = 1, gto_basis_set%nset
            maxso = MAX(maxso, gto_basis_set%npgf(iset)* &
                        nsoset(gto_basis_set%lmax(iset)))
         END DO
      END IF

      IF (PRESENT(nco_sum)) THEN
         nco_sum = 0
         DO iset = 1, gto_basis_set%nset
            nco_sum = nco_sum + gto_basis_set%npgf(iset)* &
                      ncoset(gto_basis_set%lmax(iset))
         END DO
      END IF
      IF (PRESENT(npgf_sum)) npgf_sum = SUM(gto_basis_set%npgf)
      IF (PRESENT(nshell_sum)) nshell_sum = SUM(gto_basis_set%nshell)
      IF (PRESENT(npgf_seg_sum)) THEN
         npgf_seg_sum = 0
         DO iset = 1, gto_basis_set%nset
            npgf_seg_sum = npgf_seg_sum + gto_basis_set%npgf(iset)*gto_basis_set%nshell(iset)
         END DO
      END IF

   END SUBROUTINE get_gto_basis_set

! **************************************************************************************************
!> \brief ...
!> \param gto_basis_set ...
! **************************************************************************************************
   SUBROUTINE init_aux_basis_set(gto_basis_set)

      ! Initialise a Gaussian-type orbital (GTO) basis set data set.

      ! - Creation (06.12.2000,MK)

      TYPE(gto_basis_set_type), POINTER                  :: gto_basis_set

      CHARACTER(len=*), PARAMETER :: routineN = 'init_aux_basis_set'

      INTEGER                                            :: handle

! -------------------------------------------------------------------------

      IF (.NOT. ASSOCIATED(gto_basis_set)) RETURN

      CALL timeset(routineN, handle)

      SELECT CASE (gto_basis_set%norm_type)
      CASE (0)
         ! No normalisation requested
      CASE (1)
         CALL init_norm_cgf_aux_2(gto_basis_set)
      CASE (2)
         ! WARNING this was never tested
         CALL init_norm_cgf_aux(gto_basis_set)
      CASE DEFAULT
         CPABORT("Normalization method not specified")
      END SELECT

      ! Initialise the transformation matrices "pgf" -> "cgf"
      CALL init_cphi_and_sphi(gto_basis_set)

      CALL timestop(handle)

   END SUBROUTINE init_aux_basis_set

! **************************************************************************************************
!> \brief ...
!> \param gto_basis_set ...
! **************************************************************************************************
   SUBROUTINE init_cphi_and_sphi(gto_basis_set)

      ! Initialise the matrices for the transformation of primitive Cartesian
      ! Gaussian-type functions to contracted Cartesian (cphi) and spherical
      ! (sphi) Gaussian-type functions.

      ! - Creation (20.09.2000,MK)

      TYPE(gto_basis_set_type), INTENT(INOUT)            :: gto_basis_set

      CHARACTER(len=*), PARAMETER :: routineN = 'init_cphi_and_sphi'

      INTEGER                                            :: first_cgf, first_sgf, handle, icgf, ico, &
                                                            ipgf, iset, ishell, l, last_sgf, lmax, &
                                                            lmin, n, n1, n2, ncgf, nn, nn1, nn2, &
                                                            npgf, nsgf

! -------------------------------------------------------------------------
! Build the Cartesian transformation matrix "cphi"

      CALL timeset(routineN, handle)

      gto_basis_set%cphi = 0.0_dp
      DO iset = 1, gto_basis_set%nset
         n = ncoset(gto_basis_set%lmax(iset))
         DO ishell = 1, gto_basis_set%nshell(iset)
            DO icgf = gto_basis_set%first_cgf(ishell, iset), &
               gto_basis_set%last_cgf(ishell, iset)
               ico = coset(gto_basis_set%lx(icgf), &
                           gto_basis_set%ly(icgf), &
                           gto_basis_set%lz(icgf))
               DO ipgf = 1, gto_basis_set%npgf(iset)
                  gto_basis_set%cphi(ico, icgf) = gto_basis_set%norm_cgf(icgf)* &
                                                  gto_basis_set%gcc(ipgf, ishell, iset)
                  ico = ico + n
               END DO
            END DO
         END DO
      END DO

      ! Build the spherical transformation matrix "sphi"

      n = SIZE(gto_basis_set%cphi, 1)

      gto_basis_set%sphi = 0.0_dp
      IF (n > 0) THEN
         lmax = -1
         ! Ensure proper setup of orbtramat
         DO iset = 1, gto_basis_set%nset
            DO ishell = 1, gto_basis_set%nshell(iset)
               lmax = MAX(lmax, gto_basis_set%l(ishell, iset))
            END DO
         END DO
         CALL init_spherical_harmonics(lmax, -1)

         DO iset = 1, gto_basis_set%nset
            DO ishell = 1, gto_basis_set%nshell(iset)
               l = gto_basis_set%l(ishell, iset)
               first_cgf = gto_basis_set%first_cgf(ishell, iset)
               first_sgf = gto_basis_set%first_sgf(ishell, iset)
               ncgf = nco(l)
               nsgf = nso(l)
               CALL dgemm("N", "T", n, nsgf, ncgf, &
                          1.0_dp, gto_basis_set%cphi(1, first_cgf), n, &
                          orbtramat(l)%c2s(1, 1), nsgf, &
                          0.0_dp, gto_basis_set%sphi(1, first_sgf), n)
            END DO
         END DO
      END IF

      ! Build the reduced transformation matrix "scon"
      ! This matrix transforms from cartesian primitifs to spherical contracted functions
      ! "scon" only includes primitifs (lmin -> lmax), whereas sphi is (0 -> lmax)

      n = SIZE(gto_basis_set%scon, 1)

      gto_basis_set%scon = 0.0_dp
      IF (n > 0) THEN
         DO iset = 1, gto_basis_set%nset
            lmin = gto_basis_set%lmin(iset)
            lmax = gto_basis_set%lmax(iset)
            npgf = gto_basis_set%npgf(iset)
            nn = ncoset(lmax) - ncoset(lmin - 1)
            DO ishell = 1, gto_basis_set%nshell(iset)
               first_sgf = gto_basis_set%first_sgf(ishell, iset)
               last_sgf = gto_basis_set%last_sgf(ishell, iset)
               DO ipgf = 1, npgf
                  nn1 = (ipgf - 1)*ncoset(lmax) + ncoset(lmin - 1) + 1
                  nn2 = ipgf*ncoset(lmax)
                  n1 = (ipgf - 1)*nn + 1
                  n2 = ipgf*nn
                  gto_basis_set%scon(n1:n2, first_sgf:last_sgf) = gto_basis_set%sphi(nn1:nn2, first_sgf:last_sgf)
               END DO
            END DO
         END DO
      END IF

      CALL timestop(handle)

   END SUBROUTINE init_cphi_and_sphi

! **************************************************************************************************
!> \brief ...
!> \param gto_basis_set ...
! **************************************************************************************************
   SUBROUTINE init_norm_cgf_aux(gto_basis_set)

      ! Initialise the normalization factors of the contracted Cartesian Gaussian
      ! functions, if the Gaussian functions represent charge distributions.

      ! - Creation (07.12.2000,MK)

      TYPE(gto_basis_set_type), INTENT(INOUT)            :: gto_basis_set

      INTEGER                                            :: icgf, ico, ipgf, iset, ishell, jco, &
                                                            jpgf, ll, lmax, lmin, lx, ly, lz, n, &
                                                            npgfa
      REAL(KIND=dp)                                      :: fnorm, gcca, gccb
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: ff
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: gaa
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: vv
      REAL(KIND=dp), DIMENSION(:), POINTER               :: rpgfa, zeta

! -------------------------------------------------------------------------

      n = 0
      ll = 0
      DO iset = 1, gto_basis_set%nset
         n = MAX(n, gto_basis_set%npgf(iset)*ncoset(gto_basis_set%lmax(iset)))
         ll = MAX(ll, gto_basis_set%lmax(iset))
      END DO

      ALLOCATE (gaa(n, n))
      ALLOCATE (vv(ncoset(ll), ncoset(ll), ll + ll + 1))
      ALLOCATE (ff(0:ll + ll))

      DO iset = 1, gto_basis_set%nset
         lmax = gto_basis_set%lmax(iset)
         lmin = gto_basis_set%lmin(iset)
         n = ncoset(lmax)
         npgfa = gto_basis_set%npgf(iset)
         rpgfa => gto_basis_set%pgf_radius(1:npgfa, iset)
         zeta => gto_basis_set%zet(1:npgfa, iset)
         CALL coulomb2(lmax, npgfa, zeta, rpgfa, lmin, &
                       lmax, npgfa, zeta, rpgfa, lmin, &
                       [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp, gaa, vv, ff(0:))
         DO ishell = 1, gto_basis_set%nshell(iset)
            DO icgf = gto_basis_set%first_cgf(ishell, iset), &
               gto_basis_set%last_cgf(ishell, iset)
               lx = gto_basis_set%lx(icgf)
               ly = gto_basis_set%ly(icgf)
               lz = gto_basis_set%lz(icgf)
               ico = coset(lx, ly, lz)
               fnorm = 0.0_dp
               DO ipgf = 1, npgfa
                  gcca = gto_basis_set%gcc(ipgf, ishell, iset)
                  jco = coset(lx, ly, lz)
                  DO jpgf = 1, npgfa
                     gccb = gto_basis_set%gcc(jpgf, ishell, iset)
                     fnorm = fnorm + gcca*gccb*gaa(ico, jco)
                     jco = jco + n
                  END DO
                  ico = ico + n
               END DO
               gto_basis_set%norm_cgf(icgf) = 1.0_dp/SQRT(fnorm)
            END DO
         END DO
      END DO

      DEALLOCATE (vv, ff)
      DEALLOCATE (gaa)

   END SUBROUTINE init_norm_cgf_aux

! **************************************************************************************************
!> \brief ...
!> \param gto_basis_set ...
! **************************************************************************************************
   ELEMENTAL SUBROUTINE init_norm_cgf_aux_2(gto_basis_set)

      ! Initialise the normalization factors of the auxiliary Cartesian Gaussian
      ! functions (Kim-Gordon polarization basis) Norm = 1.

      ! - Creation (07.12.2000,GT)

      TYPE(gto_basis_set_type), INTENT(INOUT)            :: gto_basis_set

      INTEGER                                            :: icgf, iset, ishell

      DO iset = 1, gto_basis_set%nset
         DO ishell = 1, gto_basis_set%nshell(iset)
            DO icgf = gto_basis_set%first_cgf(ishell, iset), &
               gto_basis_set%last_cgf(ishell, iset)
               gto_basis_set%norm_cgf(icgf) = 1.0_dp
            END DO
         END DO
      END DO

   END SUBROUTINE init_norm_cgf_aux_2

! **************************************************************************************************
!> \brief Initialise the normalization factors of the contracted Cartesian Gaussian functions.
!> \param gto_basis_set ...
!> \author MK
! **************************************************************************************************
   ELEMENTAL SUBROUTINE init_norm_cgf_orb(gto_basis_set)

      TYPE(gto_basis_set_type), INTENT(INOUT)            :: gto_basis_set

      INTEGER                                            :: icgf, ipgf, iset, ishell, jpgf, l, lx, &
                                                            ly, lz
      REAL(KIND=dp)                                      :: expzet, fnorm, gcca, gccb, prefac, zeta, &
                                                            zetb

      DO iset = 1, gto_basis_set%nset
         DO ishell = 1, gto_basis_set%nshell(iset)

            l = gto_basis_set%l(ishell, iset)

            expzet = 0.5_dp*REAL(2*l + 3, dp)

            fnorm = 0.0_dp

            DO ipgf = 1, gto_basis_set%npgf(iset)
               gcca = gto_basis_set%gcc(ipgf, ishell, iset)
               zeta = gto_basis_set%zet(ipgf, iset)
               DO jpgf = 1, gto_basis_set%npgf(iset)
                  gccb = gto_basis_set%gcc(jpgf, ishell, iset)
                  zetb = gto_basis_set%zet(jpgf, iset)
                  fnorm = fnorm + gcca*gccb/(zeta + zetb)**expzet
               END DO
            END DO

            fnorm = 0.5_dp**l*pi**1.5_dp*fnorm

            DO icgf = gto_basis_set%first_cgf(ishell, iset), &
               gto_basis_set%last_cgf(ishell, iset)
               lx = gto_basis_set%lx(icgf)
               ly = gto_basis_set%ly(icgf)
               lz = gto_basis_set%lz(icgf)
               prefac = dfac(2*lx - 1)*dfac(2*ly - 1)*dfac(2*lz - 1)
               gto_basis_set%norm_cgf(icgf) = 1.0_dp/SQRT(prefac*fnorm)
            END DO

         END DO
      END DO

   END SUBROUTINE init_norm_cgf_orb

! **************************************************************************************************
!> \brief Initialise the normalization factors of the contracted Cartesian Gaussian
!>        functions used for frozen density representation.
!> \param gto_basis_set ...
!> \author GT
! **************************************************************************************************
   ELEMENTAL SUBROUTINE init_norm_cgf_orb_den(gto_basis_set)

      TYPE(gto_basis_set_type), INTENT(INOUT)            :: gto_basis_set

      INTEGER                                            :: icgf, ipgf, iset, ishell, l
      REAL(KIND=dp)                                      :: expzet, gcca, prefac, zeta

      DO iset = 1, gto_basis_set%nset
         DO ishell = 1, gto_basis_set%nshell(iset)
            l = gto_basis_set%l(ishell, iset)
            expzet = 0.5_dp*REAL(2*l + 3, dp)
            prefac = (1.0_dp/pi)**1.5_dp
            DO ipgf = 1, gto_basis_set%npgf(iset)
               gcca = gto_basis_set%gcc(ipgf, ishell, iset)
               zeta = gto_basis_set%zet(ipgf, iset)
               gto_basis_set%gcc(ipgf, ishell, iset) = prefac*zeta**expzet*gcca
            END DO
            DO icgf = gto_basis_set%first_cgf(ishell, iset), &
               gto_basis_set%last_cgf(ishell, iset)
               gto_basis_set%norm_cgf(icgf) = 1.0_dp
            END DO
         END DO
      END DO

   END SUBROUTINE init_norm_cgf_orb_den

! **************************************************************************************************
!> \brief Initialise a Gaussian-type orbital (GTO) basis set data set.
!> \param gto_basis_set ...
!> \author MK
! **************************************************************************************************
   SUBROUTINE init_orb_basis_set(gto_basis_set)

      TYPE(gto_basis_set_type), POINTER                  :: gto_basis_set

      CHARACTER(len=*), PARAMETER :: routineN = 'init_orb_basis_set'

      INTEGER                                            :: handle

! -------------------------------------------------------------------------

      IF (.NOT. ASSOCIATED(gto_basis_set)) RETURN

      CALL timeset(routineN, handle)

      SELECT CASE (gto_basis_set%norm_type)
      CASE (0)
         ! No normalisation requested
      CASE (1)
         CALL init_norm_cgf_orb_den(gto_basis_set)
      CASE (2)
         ! Normalise the primitive Gaussian functions
         CALL normalise_gcc_orb(gto_basis_set)
         ! Compute the normalization factors of the contracted Gaussian-type
         ! functions
         CALL init_norm_cgf_orb(gto_basis_set)
      CASE (3)
         CALL init_norm_cgf_orb(gto_basis_set)
      CASE DEFAULT
         CPABORT("Normalization method not specified")
      END SELECT

      ! Initialise the transformation matrices "pgf" -> "cgf"

      CALL init_cphi_and_sphi(gto_basis_set)

      CALL timestop(handle)

   END SUBROUTINE init_orb_basis_set

! **************************************************************************************************
!> \brief Normalise the primitive Cartesian Gaussian functions. The normalization
!>        factor is included in the Gaussian contraction coefficients.
!> \param gto_basis_set ...
!> \author MK
! **************************************************************************************************
   SUBROUTINE normalise_gcc_orb(gto_basis_set)

      TYPE(gto_basis_set_type), POINTER                  :: gto_basis_set

      INTEGER                                            :: ipgf, iset, ishell, l
      REAL(KIND=dp)                                      :: expzet, gcca, prefac, zeta

      DO iset = 1, gto_basis_set%nset
         DO ishell = 1, gto_basis_set%nshell(iset)
            l = gto_basis_set%l(ishell, iset)
            expzet = 0.25_dp*REAL(2*l + 3, dp)
            prefac = 2.0_dp**l*(2.0_dp/pi)**0.75_dp
            DO ipgf = 1, gto_basis_set%npgf(iset)
               gcca = gto_basis_set%gcc(ipgf, ishell, iset)
               zeta = gto_basis_set%zet(ipgf, iset)
               gto_basis_set%gcc(ipgf, ishell, iset) = prefac*zeta**expzet*gcca
            END DO
         END DO
      END DO

   END SUBROUTINE normalise_gcc_orb

! **************************************************************************************************
!> \brief Read a Gaussian-type orbital (GTO) basis set from the database file.
!> \param element_symbol ...
!> \param basis_set_name ...
!> \param gto_basis_set ...
!> \param para_env ...
!> \param dft_section ...
!> \author MK
! **************************************************************************************************
   SUBROUTINE read_gto_basis_set1(element_symbol, basis_set_name, gto_basis_set, &
                                  para_env, dft_section)

      CHARACTER(LEN=*), INTENT(IN)                       :: element_symbol, basis_set_name
      TYPE(gto_basis_set_type), INTENT(INOUT)            :: gto_basis_set
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(section_vals_type), POINTER                   :: dft_section

      CHARACTER(LEN=240)                                 :: line
      CHARACTER(LEN=242)                                 :: line2
      CHARACTER(len=default_path_length)                 :: basis_set_file_name, tmp
      CHARACTER(LEN=default_path_length), DIMENSION(:), &
         POINTER                                         :: cbasis
      CHARACTER(LEN=LEN(basis_set_name))                 :: bsname
      CHARACTER(LEN=LEN(basis_set_name)+2)               :: bsname2
      CHARACTER(LEN=LEN(element_symbol))                 :: symbol
      CHARACTER(LEN=LEN(element_symbol)+2)               :: symbol2
      INTEGER :: i, ibasis, ico, ipgf, irep, iset, ishell, lshell, m, maxco, maxl, maxpgf, &
         maxshell, nbasis, ncgf, nmin, nset, nsgf, sort_method, strlen1, strlen2
      INTEGER, DIMENSION(:), POINTER                     :: lmax, lmin, npgf, nshell
      INTEGER, DIMENSION(:, :), POINTER                  :: l, n
      LOGICAL                                            :: basis_found, found, match
      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zet
      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: gcc
      TYPE(cp_parser_type)                               :: parser

      line = ""
      line2 = ""
      symbol = ""
      symbol2 = ""
      bsname = ""
      bsname2 = ""

      nbasis = 1

      gto_basis_set%name = basis_set_name
      gto_basis_set%aliases = basis_set_name
      CALL section_vals_val_get(dft_section, "BASIS_SET_FILE_NAME", &
                                n_rep_val=nbasis)
      ALLOCATE (cbasis(nbasis))
      DO ibasis = 1, nbasis
         CALL section_vals_val_get(dft_section, "BASIS_SET_FILE_NAME", &
                                   i_rep_val=ibasis, c_val=cbasis(ibasis))
         basis_set_file_name = cbasis(ibasis)
         tmp = basis_set_file_name
         CALL uppercase(tmp)
         IF (INDEX(tmp, "MOLOPT") /= 0) CALL cite_reference(VandeVondele2007)
      END DO

      ! Search for the requested basis set in the basis set file
      ! until the basis set is found or the end of file is reached

      basis_found = .FALSE.
      basis_loop: DO ibasis = 1, nbasis
         IF (basis_found) EXIT basis_loop
         basis_set_file_name = cbasis(ibasis)
         CALL parser_create(parser, basis_set_file_name, para_env=para_env)

         bsname = basis_set_name
         symbol = element_symbol
         irep = 0

         tmp = basis_set_name
         CALL uppercase(tmp)
         IF (INDEX(tmp, "MOLOPT") /= 0) CALL cite_reference(VandeVondele2007)

         nset = 0
         maxshell = 0
         maxpgf = 0
         maxco = 0
         ncgf = 0
         nsgf = 0
         gto_basis_set%nset = nset
         gto_basis_set%ncgf = ncgf
         gto_basis_set%nsgf = nsgf
         CALL reallocate(gto_basis_set%lmax, 1, nset)
         CALL reallocate(gto_basis_set%lmin, 1, nset)
         CALL reallocate(gto_basis_set%npgf, 1, nset)
         CALL reallocate(gto_basis_set%nshell, 1, nset)
         CALL reallocate(gto_basis_set%n, 1, maxshell, 1, nset)
         CALL reallocate(gto_basis_set%l, 1, maxshell, 1, nset)
         CALL reallocate(gto_basis_set%zet, 1, maxpgf, 1, nset)
         CALL reallocate(gto_basis_set%gcc, 1, maxpgf, 1, maxshell, 1, nset)
         CALL reallocate(gto_basis_set%set_radius, 1, nset)
         CALL reallocate(gto_basis_set%pgf_radius, 1, maxpgf, 1, nset)
         CALL reallocate(gto_basis_set%first_cgf, 1, maxshell, 1, nset)
         CALL reallocate(gto_basis_set%first_sgf, 1, maxshell, 1, nset)
         CALL reallocate(gto_basis_set%last_cgf, 1, maxshell, 1, nset)
         CALL reallocate(gto_basis_set%last_sgf, 1, maxshell, 1, nset)
         CALL reallocate(gto_basis_set%ncgf_set, 1, nset)
         CALL reallocate(gto_basis_set%nsgf_set, 1, nset)
         CALL reallocate(gto_basis_set%cphi, 1, maxco, 1, ncgf)
         CALL reallocate(gto_basis_set%sphi, 1, maxco, 1, nsgf)
         CALL reallocate(gto_basis_set%scon, 1, maxco, 1, nsgf)
         CALL reallocate(gto_basis_set%lx, 1, ncgf)
         CALL reallocate(gto_basis_set%ly, 1, ncgf)
         CALL reallocate(gto_basis_set%lz, 1, ncgf)
         CALL reallocate(gto_basis_set%m, 1, nsgf)
         CALL reallocate(gto_basis_set%norm_cgf, 1, ncgf)

         IF (tmp /= "NONE") THEN
            search_loop: DO

               CALL parser_search_string(parser, TRIM(bsname), .TRUE., found, line)
               IF (found) THEN
                  CALL uppercase(symbol)
                  CALL uppercase(bsname)

                  match = .FALSE.
                  CALL uppercase(line)
                  ! Check both the element symbol and the basis set name
                  line2 = " "//line//" "
                  symbol2 = " "//TRIM(symbol)//" "
                  bsname2 = " "//TRIM(bsname)//" "
                  strlen1 = LEN_TRIM(symbol2) + 1
                  strlen2 = LEN_TRIM(bsname2) + 1

                  IF ((INDEX(line2, symbol2(:strlen1)) > 0) .AND. &
                      (INDEX(line2, bsname2(:strlen2)) > 0)) match = .TRUE.
                  IF (match) THEN
                     ! copy all names into aliases field
                     i = INDEX(line2, symbol2(:strlen1))
                     i = i + 1 + INDEX(line2(i + 1:), " ")
                     gto_basis_set%aliases = line2(i:)

                     NULLIFY (gcc, l, lmax, lmin, n, npgf, nshell, zet)
                     ! Read the basis set information
                     CALL parser_get_object(parser, nset, newline=.TRUE.)

                     CALL reallocate(npgf, 1, nset)
                     CALL reallocate(nshell, 1, nset)
                     CALL reallocate(lmax, 1, nset)
                     CALL reallocate(lmin, 1, nset)
                     CALL reallocate(n, 1, 1, 1, nset)

                     maxl = 0
                     maxpgf = 0
                     maxshell = 0

                     DO iset = 1, nset
                        CALL parser_get_object(parser, n(1, iset), newline=.TRUE.)
                        CALL parser_get_object(parser, lmin(iset))
                        CALL parser_get_object(parser, lmax(iset))
                        CALL parser_get_object(parser, npgf(iset))
                        maxl = MAX(maxl, lmax(iset))
                        IF (npgf(iset) > maxpgf) THEN
                           maxpgf = npgf(iset)
                           CALL reallocate(zet, 1, maxpgf, 1, nset)
                           CALL reallocate(gcc, 1, maxpgf, 1, maxshell, 1, nset)
                        END IF
                        nshell(iset) = 0
                        DO lshell = lmin(iset), lmax(iset)
                           nmin = n(1, iset) + lshell - lmin(iset)
                           CALL parser_get_object(parser, ishell)
                           nshell(iset) = nshell(iset) + ishell
                           IF (nshell(iset) > maxshell) THEN
                              maxshell = nshell(iset)
                              CALL reallocate(n, 1, maxshell, 1, nset)
                              CALL reallocate(l, 1, maxshell, 1, nset)
                              CALL reallocate(gcc, 1, maxpgf, 1, maxshell, 1, nset)
                           END IF
                           DO i = 1, ishell
                              n(nshell(iset) - ishell + i, iset) = nmin + i - 1
                              l(nshell(iset) - ishell + i, iset) = lshell
                           END DO
                        END DO
                        DO ipgf = 1, npgf(iset)
                           CALL parser_get_object(parser, zet(ipgf, iset), newline=.TRUE.)
                           DO ishell = 1, nshell(iset)
                              CALL parser_get_object(parser, gcc(ipgf, ishell, iset))
                           END DO
                        END DO
                     END DO

                     ! Maximum angular momentum quantum number of the atomic kind

                     CALL init_orbital_pointers(maxl)

                     ! Allocate the global variables

                     gto_basis_set%nset = nset
                     CALL reallocate(gto_basis_set%lmax, 1, nset)
                     CALL reallocate(gto_basis_set%lmin, 1, nset)
                     CALL reallocate(gto_basis_set%npgf, 1, nset)
                     CALL reallocate(gto_basis_set%nshell, 1, nset)
                     CALL reallocate(gto_basis_set%n, 1, maxshell, 1, nset)
                     CALL reallocate(gto_basis_set%l, 1, maxshell, 1, nset)
                     CALL reallocate(gto_basis_set%zet, 1, maxpgf, 1, nset)
                     CALL reallocate(gto_basis_set%gcc, 1, maxpgf, 1, maxshell, 1, nset)

                     ! Copy the basis set information into the data structure

                     DO iset = 1, nset
                        gto_basis_set%lmax(iset) = lmax(iset)
                        gto_basis_set%lmin(iset) = lmin(iset)
                        gto_basis_set%npgf(iset) = npgf(iset)
                        gto_basis_set%nshell(iset) = nshell(iset)
                        DO ishell = 1, nshell(iset)
                           gto_basis_set%n(ishell, iset) = n(ishell, iset)
                           gto_basis_set%l(ishell, iset) = l(ishell, iset)
                           DO ipgf = 1, npgf(iset)
                              gto_basis_set%gcc(ipgf, ishell, iset) = gcc(ipgf, ishell, iset)
                           END DO
                        END DO
                        DO ipgf = 1, npgf(iset)
                           gto_basis_set%zet(ipgf, iset) = zet(ipgf, iset)
                        END DO
                     END DO

                     ! Initialise the depending atomic kind information

                     CALL reallocate(gto_basis_set%set_radius, 1, nset)
                     CALL reallocate(gto_basis_set%pgf_radius, 1, maxpgf, 1, nset)
                     CALL reallocate(gto_basis_set%first_cgf, 1, maxshell, 1, nset)
                     CALL reallocate(gto_basis_set%first_sgf, 1, maxshell, 1, nset)
                     CALL reallocate(gto_basis_set%last_cgf, 1, maxshell, 1, nset)
                     CALL reallocate(gto_basis_set%last_sgf, 1, maxshell, 1, nset)
                     CALL reallocate(gto_basis_set%ncgf_set, 1, nset)
                     CALL reallocate(gto_basis_set%nsgf_set, 1, nset)

                     maxco = 0
                     ncgf = 0
                     nsgf = 0

                     DO iset = 1, nset
                        gto_basis_set%ncgf_set(iset) = 0
                        gto_basis_set%nsgf_set(iset) = 0
                        DO ishell = 1, nshell(iset)
                           lshell = gto_basis_set%l(ishell, iset)
                           gto_basis_set%first_cgf(ishell, iset) = ncgf + 1
                           ncgf = ncgf + nco(lshell)
                           gto_basis_set%last_cgf(ishell, iset) = ncgf
                           gto_basis_set%ncgf_set(iset) = &
                              gto_basis_set%ncgf_set(iset) + nco(lshell)
                           gto_basis_set%first_sgf(ishell, iset) = nsgf + 1
                           nsgf = nsgf + nso(lshell)
                           gto_basis_set%last_sgf(ishell, iset) = nsgf
                           gto_basis_set%nsgf_set(iset) = &
                              gto_basis_set%nsgf_set(iset) + nso(lshell)
                        END DO
                        maxco = MAX(maxco, npgf(iset)*ncoset(lmax(iset)))
                     END DO

                     gto_basis_set%ncgf = ncgf
                     gto_basis_set%nsgf = nsgf

                     CALL reallocate(gto_basis_set%cphi, 1, maxco, 1, ncgf)
                     CALL reallocate(gto_basis_set%sphi, 1, maxco, 1, nsgf)
                     CALL reallocate(gto_basis_set%scon, 1, maxco, 1, nsgf)
                     CALL reallocate(gto_basis_set%lx, 1, ncgf)
                     CALL reallocate(gto_basis_set%ly, 1, ncgf)
                     CALL reallocate(gto_basis_set%lz, 1, ncgf)
                     CALL reallocate(gto_basis_set%m, 1, nsgf)
                     CALL reallocate(gto_basis_set%norm_cgf, 1, ncgf)
                     ALLOCATE (gto_basis_set%cgf_symbol(ncgf))

                     ALLOCATE (gto_basis_set%sgf_symbol(nsgf))

                     ncgf = 0
                     nsgf = 0

                     DO iset = 1, nset
                        DO ishell = 1, nshell(iset)
                           lshell = gto_basis_set%l(ishell, iset)
                           DO ico = ncoset(lshell - 1) + 1, ncoset(lshell)
                              ncgf = ncgf + 1
                              gto_basis_set%lx(ncgf) = indco(1, ico)
                              gto_basis_set%ly(ncgf) = indco(2, ico)
                              gto_basis_set%lz(ncgf) = indco(3, ico)
                              gto_basis_set%cgf_symbol(ncgf) = &
                                 cgf_symbol(n(ishell, iset), [gto_basis_set%lx(ncgf), &
                                                              gto_basis_set%ly(ncgf), &
                                                              gto_basis_set%lz(ncgf)])
                           END DO
                           DO m = -lshell, lshell
                              nsgf = nsgf + 1
                              gto_basis_set%m(nsgf) = m
                              gto_basis_set%sgf_symbol(nsgf) = &
                                 sgf_symbol(n(ishell, iset), lshell, m)
                           END DO
                        END DO
                     END DO

                     DEALLOCATE (gcc, l, lmax, lmin, n, npgf, nshell, zet)

                     basis_found = .TRUE.
                     EXIT search_loop
                  END IF
               ELSE
                  EXIT search_loop
               END IF
            END DO search_loop
         ELSE
            match = .FALSE.
            ALLOCATE (gto_basis_set%cgf_symbol(ncgf))
            ALLOCATE (gto_basis_set%sgf_symbol(nsgf))
         END IF

         CALL parser_release(parser)

      END DO basis_loop

      IF (tmp /= "NONE") THEN
         IF (.NOT. basis_found) THEN
            basis_set_file_name = ""
            DO ibasis = 1, nbasis
               basis_set_file_name = TRIM(basis_set_file_name)//"<"//TRIM(cbasis(ibasis))//"> "
            END DO
            CALL cp_abort(__LOCATION__, &
                          "The requested basis set <"//TRIM(bsname)// &
                          "> for element <"//TRIM(symbol)//"> was not "// &
                          "found in the basis set files "// &
                          TRIM(basis_set_file_name))
         END IF
      END IF
      DEALLOCATE (cbasis)

      CALL section_vals_val_get(dft_section, "SORT_BASIS", i_val=sort_method)
      CALL sort_gto_basis_set(gto_basis_set, sort_method)

   END SUBROUTINE read_gto_basis_set1

! **************************************************************************************************
!> \brief Read a Gaussian-type orbital (GTO) basis set from the database file.
!> \param element_symbol ...
!> \param basis_type ...
!> \param gto_basis_set ...
!> \param basis_section ...
!> \param irep ...
!> \param dft_section ...
!> \author MK
! **************************************************************************************************
   SUBROUTINE read_gto_basis_set2(element_symbol, basis_type, gto_basis_set, &
                                  basis_section, irep, dft_section)

      CHARACTER(LEN=*), INTENT(IN)                       :: element_symbol
      CHARACTER(LEN=*), INTENT(INOUT)                    :: basis_type
      TYPE(gto_basis_set_type), INTENT(INOUT)            :: gto_basis_set
      TYPE(section_vals_type), OPTIONAL, POINTER         :: basis_section
      INTEGER                                            :: irep
      TYPE(section_vals_type), OPTIONAL, POINTER         :: dft_section

      CHARACTER(len=20*default_string_length)            :: line_att
      CHARACTER(LEN=240)                                 :: line
      CHARACTER(LEN=242)                                 :: line2
      CHARACTER(LEN=default_path_length)                 :: bsname, bsname2
      CHARACTER(LEN=LEN(element_symbol))                 :: symbol
      CHARACTER(LEN=LEN(element_symbol)+2)               :: symbol2
      INTEGER                                            :: i, ico, ipgf, iset, ishell, lshell, m, &
                                                            maxco, maxl, maxpgf, maxshell, ncgf, &
                                                            nmin, nset, nsgf, sort_method
      INTEGER, DIMENSION(:), POINTER                     :: lmax, lmin, npgf, nshell
      INTEGER, DIMENSION(:, :), POINTER                  :: l, n
      LOGICAL                                            :: is_ok
      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zet
      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: gcc
      TYPE(cp_sll_val_type), POINTER                     :: list
      TYPE(val_type), POINTER                            :: val

      line = ""
      line2 = ""
      symbol = ""
      symbol2 = ""
      bsname = ""
      bsname2 = ""

      gto_basis_set%name = " "
      gto_basis_set%aliases = " "

      bsname = " "
      symbol = element_symbol

      nset = 0
      maxshell = 0
      maxpgf = 0
      maxco = 0
      ncgf = 0
      nsgf = 0
      gto_basis_set%nset = nset
      gto_basis_set%ncgf = ncgf
      gto_basis_set%nsgf = nsgf
      CALL reallocate(gto_basis_set%lmax, 1, nset)
      CALL reallocate(gto_basis_set%lmin, 1, nset)
      CALL reallocate(gto_basis_set%npgf, 1, nset)
      CALL reallocate(gto_basis_set%nshell, 1, nset)
      CALL reallocate(gto_basis_set%n, 1, maxshell, 1, nset)
      CALL reallocate(gto_basis_set%l, 1, maxshell, 1, nset)
      CALL reallocate(gto_basis_set%zet, 1, maxpgf, 1, nset)
      CALL reallocate(gto_basis_set%gcc, 1, maxpgf, 1, maxshell, 1, nset)
      CALL reallocate(gto_basis_set%set_radius, 1, nset)
      CALL reallocate(gto_basis_set%pgf_radius, 1, maxpgf, 1, nset)
      CALL reallocate(gto_basis_set%first_cgf, 1, maxshell, 1, nset)
      CALL reallocate(gto_basis_set%first_sgf, 1, maxshell, 1, nset)
      CALL reallocate(gto_basis_set%last_cgf, 1, maxshell, 1, nset)
      CALL reallocate(gto_basis_set%last_sgf, 1, maxshell, 1, nset)
      CALL reallocate(gto_basis_set%ncgf_set, 1, nset)
      CALL reallocate(gto_basis_set%nsgf_set, 1, nset)
      CALL reallocate(gto_basis_set%cphi, 1, maxco, 1, ncgf)
      CALL reallocate(gto_basis_set%sphi, 1, maxco, 1, nsgf)
      CALL reallocate(gto_basis_set%scon, 1, maxco, 1, nsgf)
      CALL reallocate(gto_basis_set%lx, 1, ncgf)
      CALL reallocate(gto_basis_set%ly, 1, ncgf)
      CALL reallocate(gto_basis_set%lz, 1, ncgf)
      CALL reallocate(gto_basis_set%m, 1, nsgf)
      CALL reallocate(gto_basis_set%norm_cgf, 1, ncgf)

      basis_type = ""
      CALL section_vals_val_get(basis_section, "_SECTION_PARAMETERS_", i_rep_section=irep, c_val=basis_type)
      IF (basis_type == "Orbital") basis_type = "ORB"

      NULLIFY (list, val)
      CALL section_vals_list_get(basis_section, "_DEFAULT_KEYWORD_", i_rep_section=irep, list=list)
      CALL uppercase(symbol)
      CALL uppercase(bsname)

      NULLIFY (gcc, l, lmax, lmin, n, npgf, nshell, zet)
      ! Read the basis set information
      is_ok = cp_sll_val_next(list, val)
      IF (.NOT. is_ok) CPABORT("Error reading the Basis set from input file!")
      CALL val_get(val, c_val=line_att)
      READ (line_att, *) nset

      CALL reallocate(npgf, 1, nset)
      CALL reallocate(nshell, 1, nset)
      CALL reallocate(lmax, 1, nset)
      CALL reallocate(lmin, 1, nset)
      CALL reallocate(n, 1, 1, 1, nset)

      maxl = 0
      maxpgf = 0
      maxshell = 0

      DO iset = 1, nset
         is_ok = cp_sll_val_next(list, val)
         IF (.NOT. is_ok) CPABORT("Error reading the Basis set from input file!")
         CALL val_get(val, c_val=line_att)
         READ (line_att, *) n(1, iset)
         CALL remove_word(line_att)
         READ (line_att, *) lmin(iset)
         CALL remove_word(line_att)
         READ (line_att, *) lmax(iset)
         CALL remove_word(line_att)
         READ (line_att, *) npgf(iset)
         CALL remove_word(line_att)
         maxl = MAX(maxl, lmax(iset))
         IF (npgf(iset) > maxpgf) THEN
            maxpgf = npgf(iset)
            CALL reallocate(zet, 1, maxpgf, 1, nset)
            CALL reallocate(gcc, 1, maxpgf, 1, maxshell, 1, nset)
         END IF
         nshell(iset) = 0
         DO lshell = lmin(iset), lmax(iset)
            nmin = n(1, iset) + lshell - lmin(iset)
            READ (line_att, *) ishell
            CALL remove_word(line_att)
            nshell(iset) = nshell(iset) + ishell
            IF (nshell(iset) > maxshell) THEN
               maxshell = nshell(iset)
               CALL reallocate(n, 1, maxshell, 1, nset)
               CALL reallocate(l, 1, maxshell, 1, nset)
               CALL reallocate(gcc, 1, maxpgf, 1, maxshell, 1, nset)
            END IF
            DO i = 1, ishell
               n(nshell(iset) - ishell + i, iset) = nmin + i - 1
               l(nshell(iset) - ishell + i, iset) = lshell
            END DO
         END DO
         IF (LEN_TRIM(line_att) /= 0) &
            CPABORT("Error reading the Basis from input file!")
         DO ipgf = 1, npgf(iset)
            is_ok = cp_sll_val_next(list, val)
            IF (.NOT. is_ok) CPABORT("Error reading the Basis set from input file!")
            CALL val_get(val, c_val=line_att)
            READ (line_att, *) zet(ipgf, iset), (gcc(ipgf, ishell, iset), ishell=1, nshell(iset))
         END DO
      END DO

      ! Maximum angular momentum quantum number of the atomic kind

      CALL init_orbital_pointers(maxl)

      ! Allocate the global variables

      gto_basis_set%nset = nset
      CALL reallocate(gto_basis_set%lmax, 1, nset)
      CALL reallocate(gto_basis_set%lmin, 1, nset)
      CALL reallocate(gto_basis_set%npgf, 1, nset)
      CALL reallocate(gto_basis_set%nshell, 1, nset)
      CALL reallocate(gto_basis_set%n, 1, maxshell, 1, nset)
      CALL reallocate(gto_basis_set%l, 1, maxshell, 1, nset)
      CALL reallocate(gto_basis_set%zet, 1, maxpgf, 1, nset)
      CALL reallocate(gto_basis_set%gcc, 1, maxpgf, 1, maxshell, 1, nset)

      ! Copy the basis set information into the data structure

      DO iset = 1, nset
         gto_basis_set%lmax(iset) = lmax(iset)
         gto_basis_set%lmin(iset) = lmin(iset)
         gto_basis_set%npgf(iset) = npgf(iset)
         gto_basis_set%nshell(iset) = nshell(iset)
         DO ishell = 1, nshell(iset)
            gto_basis_set%n(ishell, iset) = n(ishell, iset)
            gto_basis_set%l(ishell, iset) = l(ishell, iset)
            DO ipgf = 1, npgf(iset)
               gto_basis_set%gcc(ipgf, ishell, iset) = gcc(ipgf, ishell, iset)
            END DO
         END DO
         DO ipgf = 1, npgf(iset)
            gto_basis_set%zet(ipgf, iset) = zet(ipgf, iset)
         END DO
      END DO

      ! Initialise the depending atomic kind information

      CALL reallocate(gto_basis_set%set_radius, 1, nset)
      CALL reallocate(gto_basis_set%pgf_radius, 1, maxpgf, 1, nset)
      CALL reallocate(gto_basis_set%first_cgf, 1, maxshell, 1, nset)
      CALL reallocate(gto_basis_set%first_sgf, 1, maxshell, 1, nset)
      CALL reallocate(gto_basis_set%last_cgf, 1, maxshell, 1, nset)
      CALL reallocate(gto_basis_set%last_sgf, 1, maxshell, 1, nset)
      CALL reallocate(gto_basis_set%ncgf_set, 1, nset)
      CALL reallocate(gto_basis_set%nsgf_set, 1, nset)

      maxco = 0
      ncgf = 0
      nsgf = 0

      DO iset = 1, nset
         gto_basis_set%ncgf_set(iset) = 0
         gto_basis_set%nsgf_set(iset) = 0
         DO ishell = 1, nshell(iset)
            lshell = gto_basis_set%l(ishell, iset)
            gto_basis_set%first_cgf(ishell, iset) = ncgf + 1
            ncgf = ncgf + nco(lshell)
            gto_basis_set%last_cgf(ishell, iset) = ncgf
            gto_basis_set%ncgf_set(iset) = &
               gto_basis_set%ncgf_set(iset) + nco(lshell)
            gto_basis_set%first_sgf(ishell, iset) = nsgf + 1
            nsgf = nsgf + nso(lshell)
            gto_basis_set%last_sgf(ishell, iset) = nsgf
            gto_basis_set%nsgf_set(iset) = &
               gto_basis_set%nsgf_set(iset) + nso(lshell)
         END DO
         maxco = MAX(maxco, npgf(iset)*ncoset(lmax(iset)))
      END DO

      gto_basis_set%ncgf = ncgf
      gto_basis_set%nsgf = nsgf

      CALL reallocate(gto_basis_set%cphi, 1, maxco, 1, ncgf)
      CALL reallocate(gto_basis_set%sphi, 1, maxco, 1, nsgf)
      CALL reallocate(gto_basis_set%scon, 1, maxco, 1, nsgf)
      CALL reallocate(gto_basis_set%lx, 1, ncgf)
      CALL reallocate(gto_basis_set%ly, 1, ncgf)
      CALL reallocate(gto_basis_set%lz, 1, ncgf)
      CALL reallocate(gto_basis_set%m, 1, nsgf)
      CALL reallocate(gto_basis_set%norm_cgf, 1, ncgf)
      ALLOCATE (gto_basis_set%cgf_symbol(ncgf))

      ALLOCATE (gto_basis_set%sgf_symbol(nsgf))

      ncgf = 0
      nsgf = 0

      DO iset = 1, nset
         DO ishell = 1, nshell(iset)
            lshell = gto_basis_set%l(ishell, iset)
            DO ico = ncoset(lshell - 1) + 1, ncoset(lshell)
               ncgf = ncgf + 1
               gto_basis_set%lx(ncgf) = indco(1, ico)
               gto_basis_set%ly(ncgf) = indco(2, ico)
               gto_basis_set%lz(ncgf) = indco(3, ico)
               gto_basis_set%cgf_symbol(ncgf) = &
                  cgf_symbol(n(ishell, iset), [gto_basis_set%lx(ncgf), &
                                               gto_basis_set%ly(ncgf), &
                                               gto_basis_set%lz(ncgf)])
            END DO
            DO m = -lshell, lshell
               nsgf = nsgf + 1
               gto_basis_set%m(nsgf) = m
               gto_basis_set%sgf_symbol(nsgf) = &
                  sgf_symbol(n(ishell, iset), lshell, m)
            END DO
         END DO
      END DO

      DEALLOCATE (gcc, l, lmax, lmin, n, npgf, nshell, zet)

      IF (PRESENT(dft_section)) THEN
         CALL section_vals_val_get(dft_section, "SORT_BASIS", i_val=sort_method)
         CALL sort_gto_basis_set(gto_basis_set, sort_method)
      END IF

   END SUBROUTINE read_gto_basis_set2

! **************************************************************************************************
!> \brief Set the components of Gaussian-type orbital (GTO) basis set data set.
!> \param gto_basis_set ...
!> \param name ...
!> \param aliases ...
!> \param norm_type ...
!> \param kind_radius ...
!> \param ncgf ...
!> \param nset ...
!> \param nsgf ...
!> \param cgf_symbol ...
!> \param sgf_symbol ...
!> \param norm_cgf ...
!> \param set_radius ...
!> \param lmax ...
!> \param lmin ...
!> \param lx ...
!> \param ly ...
!> \param lz ...
!> \param m ...
!> \param ncgf_set ...
!> \param npgf ...
!> \param nsgf_set ...
!> \param nshell ...
!> \param cphi ...
!> \param pgf_radius ...
!> \param sphi ...
!> \param scon ...
!> \param zet ...
!> \param first_cgf ...
!> \param first_sgf ...
!> \param l ...
!> \param last_cgf ...
!> \param last_sgf ...
!> \param n ...
!> \param gcc ...
!> \param short_kind_radius ...
!> \author MK
! **************************************************************************************************
   SUBROUTINE set_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, &
                                nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, &
                                lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, &
                                cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, &
                                last_cgf, last_sgf, n, gcc, short_kind_radius)

      TYPE(gto_basis_set_type), INTENT(INOUT)            :: gto_basis_set
      CHARACTER(LEN=default_string_length), INTENT(IN), &
         OPTIONAL                                        :: name, aliases
      INTEGER, INTENT(IN), OPTIONAL                      :: norm_type
      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: kind_radius
      INTEGER, INTENT(IN), OPTIONAL                      :: ncgf, nset, nsgf
      CHARACTER(LEN=12), DIMENSION(:), OPTIONAL, POINTER :: cgf_symbol
      CHARACTER(LEN=6), DIMENSION(:), OPTIONAL, POINTER  :: sgf_symbol
      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: norm_cgf, set_radius
      INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: lmax, lmin, lx, ly, lz, m, ncgf_set, &
                                                            npgf, nsgf_set, nshell
      REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: cphi, pgf_radius, sphi, scon, zet
      INTEGER, DIMENSION(:, :), OPTIONAL, POINTER        :: first_cgf, first_sgf, l, last_cgf, &
                                                            last_sgf, n
      REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
         POINTER                                         :: gcc
      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: short_kind_radius

      IF (PRESENT(name)) gto_basis_set%name = name
      IF (PRESENT(aliases)) gto_basis_set%aliases = aliases
      IF (PRESENT(norm_type)) gto_basis_set%norm_type = norm_type
      IF (PRESENT(kind_radius)) gto_basis_set%kind_radius = kind_radius
      IF (PRESENT(short_kind_radius)) gto_basis_set%short_kind_radius = short_kind_radius
      IF (PRESENT(ncgf)) gto_basis_set%ncgf = ncgf
      IF (PRESENT(nset)) gto_basis_set%nset = nset
      IF (PRESENT(nsgf)) gto_basis_set%nsgf = nsgf
      IF (PRESENT(cgf_symbol)) gto_basis_set%cgf_symbol(:) = cgf_symbol(:)
      IF (PRESENT(sgf_symbol)) gto_basis_set%sgf_symbol(:) = sgf_symbol(:)
      IF (PRESENT(norm_cgf)) gto_basis_set%norm_cgf(:) = norm_cgf(:)
      IF (PRESENT(set_radius)) gto_basis_set%set_radius(:) = set_radius(:)
      IF (PRESENT(lmax)) gto_basis_set%lmax(:) = lmax(:)
      IF (PRESENT(lmin)) gto_basis_set%lmin(:) = lmin(:)
      IF (PRESENT(lx)) gto_basis_set%lx(:) = lx(:)
      IF (PRESENT(ly)) gto_basis_set%ly(:) = ly(:)
      IF (PRESENT(lz)) gto_basis_set%lz(:) = lz(:)
      IF (PRESENT(m)) gto_basis_set%m(:) = m(:)
      IF (PRESENT(ncgf_set)) gto_basis_set%ncgf_set(:) = ncgf_set(:)
      IF (PRESENT(npgf)) gto_basis_set%npgf(:) = npgf(:)
      IF (PRESENT(nsgf_set)) gto_basis_set%nsgf_set(:) = nsgf_set(:)
      IF (PRESENT(nshell)) gto_basis_set%nshell(:) = nshell(:)
      IF (PRESENT(cphi)) gto_basis_set%cphi(:, :) = cphi(:, :)
      IF (PRESENT(pgf_radius)) gto_basis_set%pgf_radius(:, :) = pgf_radius(:, :)
      IF (PRESENT(sphi)) gto_basis_set%sphi(:, :) = sphi(:, :)
      IF (PRESENT(scon)) gto_basis_set%scon(:, :) = scon(:, :)
      IF (PRESENT(zet)) gto_basis_set%zet(:, :) = zet(:, :)
      IF (PRESENT(first_cgf)) gto_basis_set%first_cgf(:, :) = first_cgf(:, :)
      IF (PRESENT(first_sgf)) gto_basis_set%first_sgf(:, :) = first_sgf(:, :)
      IF (PRESENT(l)) l(:, :) = gto_basis_set%l(:, :)
      IF (PRESENT(last_cgf)) gto_basis_set%last_cgf(:, :) = last_cgf(:, :)
      IF (PRESENT(last_sgf)) gto_basis_set%last_sgf(:, :) = last_sgf(:, :)
      IF (PRESENT(n)) gto_basis_set%n(:, :) = n(:, :)
      IF (PRESENT(gcc)) gto_basis_set%gcc(:, :, :) = gcc(:, :, :)

   END SUBROUTINE set_gto_basis_set

! **************************************************************************************************
!> \brief Write a Gaussian-type orbital (GTO) basis set data set to the output unit.
!> \param gto_basis_set ...
!> \param output_unit ...
!> \param header ...
!> \author MK
! **************************************************************************************************
   SUBROUTINE write_gto_basis_set(gto_basis_set, output_unit, header)

      TYPE(gto_basis_set_type), INTENT(IN)               :: gto_basis_set
      INTEGER, INTENT(in)                                :: output_unit
      CHARACTER(len=*), OPTIONAL                         :: header

      INTEGER                                            :: ipgf, iset, ishell

      IF (output_unit > 0) THEN

         IF (PRESENT(header)) THEN
            WRITE (UNIT=output_unit, FMT="(/,T6,A,T41,A40)") &
               TRIM(header), TRIM(gto_basis_set%name)
         END IF

         WRITE (UNIT=output_unit, FMT="(/,(T8,A,T71,I10))") &
            "Number of orbital shell sets:            ", &
            gto_basis_set%nset, &
            "Number of orbital shells:                ", &
            SUM(gto_basis_set%nshell(:)), &
            "Number of primitive Cartesian functions: ", &
            SUM(gto_basis_set%npgf(:)), &
            "Number of Cartesian basis functions:     ", &
            gto_basis_set%ncgf, &
            "Number of spherical basis functions:     ", &
            gto_basis_set%nsgf, &
            "Norm type:                               ", &
            gto_basis_set%norm_type

         WRITE (UNIT=output_unit, FMT="(/,T6,A,T41,A40,/,/,T25,A)") &
            "GTO basis set information for", TRIM(gto_basis_set%name), &
            "Set   Shell     n   l            Exponent    Coefficient"

         DO iset = 1, gto_basis_set%nset
            WRITE (UNIT=output_unit, FMT="(A)") ""
            DO ishell = 1, gto_basis_set%nshell(iset)
               WRITE (UNIT=output_unit, &
                      FMT="(T25,I3,4X,I4,4X,I2,2X,I2,(T51,2F15.6))") &
                  iset, ishell, &
                  gto_basis_set%n(ishell, iset), &
                  gto_basis_set%l(ishell, iset), &
                  (gto_basis_set%zet(ipgf, iset), &
                   gto_basis_set%gcc(ipgf, ishell, iset), &
                   ipgf=1, gto_basis_set%npgf(iset))
            END DO
         END DO

      END IF

   END SUBROUTINE write_gto_basis_set

! **************************************************************************************************

! **************************************************************************************************
!> \brief Write a Gaussian-type orbital (GTO) basis set data set to the output unit.
!> \param orb_basis_set ...
!> \param output_unit ...
!> \param header ...
!> \author MK
! **************************************************************************************************
   SUBROUTINE write_orb_basis_set(orb_basis_set, output_unit, header)

      TYPE(gto_basis_set_type), INTENT(IN)               :: orb_basis_set
      INTEGER, INTENT(in)                                :: output_unit
      CHARACTER(len=*), OPTIONAL                         :: header

      INTEGER                                            :: icgf, ico, ipgf, iset, ishell

      IF (output_unit > 0) THEN
         IF (PRESENT(header)) THEN
            WRITE (UNIT=output_unit, FMT="(/,T6,A,T41,A40)") &
               TRIM(header), TRIM(orb_basis_set%name)
         END IF

         WRITE (UNIT=output_unit, FMT="(/,(T8,A,T71,I10))") &
            "Number of orbital shell sets:            ", &
            orb_basis_set%nset, &
            "Number of orbital shells:                ", &
            SUM(orb_basis_set%nshell(:)), &
            "Number of primitive Cartesian functions: ", &
            SUM(orb_basis_set%npgf(:)), &
            "Number of Cartesian basis functions:     ", &
            orb_basis_set%ncgf, &
            "Number of spherical basis functions:     ", &
            orb_basis_set%nsgf, &
            "Norm type:                               ", &
            orb_basis_set%norm_type

         WRITE (UNIT=output_unit, FMT="(/,T8,A,/,/,T25,A)") &
            "Normalised Cartesian orbitals:", &
            "Set   Shell   Orbital            Exponent    Coefficient"

         icgf = 0

         DO iset = 1, orb_basis_set%nset
            DO ishell = 1, orb_basis_set%nshell(iset)
               WRITE (UNIT=output_unit, FMT="(A)") ""
               DO ico = 1, nco(orb_basis_set%l(ishell, iset))
                  icgf = icgf + 1
                  WRITE (UNIT=output_unit, &
                         FMT="(T25,I3,4X,I4,3X,A12,(T51,2F15.6))") &
                     iset, ishell, orb_basis_set%cgf_symbol(icgf), &
                     (orb_basis_set%zet(ipgf, iset), &
                      orb_basis_set%norm_cgf(icgf)* &
                      orb_basis_set%gcc(ipgf, ishell, iset), &
                      ipgf=1, orb_basis_set%npgf(iset))
               END DO
            END DO
         END DO
      END IF

   END SUBROUTINE write_orb_basis_set

! **************************************************************************************************
!> \brief ...
!> \param gto_basis_set ...
!> \param output_unit ...
!> \author JGH
! **************************************************************************************************
   SUBROUTINE dump_gto_basis_set(gto_basis_set, output_unit)

      TYPE(gto_basis_set_type), INTENT(IN)               :: gto_basis_set
      INTEGER, INTENT(in)                                :: output_unit

      INTEGER                                            :: ipgf, iset, ishell

      IF (output_unit > 0) THEN

         WRITE (UNIT=output_unit, FMT="(/,T6,A40)") TRIM(gto_basis_set%name)
         WRITE (UNIT=output_unit, FMT="(/,T6,A40)") TRIM(gto_basis_set%aliases)
         WRITE (UNIT=output_unit, FMT="(/,T6,F12.8)") gto_basis_set%kind_radius
         WRITE (UNIT=output_unit, FMT="(/,T6,F12.8)") gto_basis_set%short_kind_radius
         WRITE (UNIT=output_unit, FMT="(/,T6,I8)") gto_basis_set%norm_type
         WRITE (UNIT=output_unit, FMT="(/,T6,3I8)") gto_basis_set%ncgf, gto_basis_set%nset, gto_basis_set%nsgf
         WRITE (UNIT=output_unit, FMT="(/,T6,6A12)") gto_basis_set%cgf_symbol
         WRITE (UNIT=output_unit, FMT="(/,T6,6A12)") gto_basis_set%sgf_symbol
         WRITE (UNIT=output_unit, FMT="(/,T6,6F12.6)") gto_basis_set%norm_cgf
         WRITE (UNIT=output_unit, FMT="(/,T6,6F12.6)") gto_basis_set%set_radius
         WRITE (UNIT=output_unit, FMT="(/,T6,12I6)") gto_basis_set%lmax
         WRITE (UNIT=output_unit, FMT="(/,T6,12I6)") gto_basis_set%lmin
         WRITE (UNIT=output_unit, FMT="(/,T6,12I6)") gto_basis_set%lx
         WRITE (UNIT=output_unit, FMT="(/,T6,12I6)") gto_basis_set%ly
         WRITE (UNIT=output_unit, FMT="(/,T6,12I6)") gto_basis_set%lz
         WRITE (UNIT=output_unit, FMT="(/,T6,12I6)") gto_basis_set%m
         WRITE (UNIT=output_unit, FMT="(/,T6,12I6)") gto_basis_set%ncgf_set
         WRITE (UNIT=output_unit, FMT="(/,T6,12I6)") gto_basis_set%nsgf_set
         WRITE (UNIT=output_unit, FMT="(/,T6,12I6)") gto_basis_set%npgf
         WRITE (UNIT=output_unit, FMT="(/,T6,12I6)") gto_basis_set%nshell

         DO iset = 1, gto_basis_set%nset
            WRITE (UNIT=output_unit, FMT="(T8,6F15.6)") &
               gto_basis_set%pgf_radius(1:gto_basis_set%npgf(iset), iset)
         END DO

         DO iset = 1, gto_basis_set%nset
            WRITE (UNIT=output_unit, FMT="(T8,6F15.6)") &
               gto_basis_set%zet(1:gto_basis_set%npgf(iset), iset)
         END DO

         DO iset = 1, gto_basis_set%nset
            DO ishell = 1, gto_basis_set%nshell(iset)
               WRITE (UNIT=output_unit, FMT="(T8,8I10)") &
                  iset, ishell, &
                  gto_basis_set%n(ishell, iset), &
                  gto_basis_set%l(ishell, iset), &
                  gto_basis_set%first_cgf(ishell, iset), &
                  gto_basis_set%last_cgf(ishell, iset), &
                  gto_basis_set%first_sgf(ishell, iset), &
                  gto_basis_set%last_sgf(ishell, iset)
            END DO
         END DO

         DO iset = 1, gto_basis_set%nset
            DO ishell = 1, gto_basis_set%nshell(iset)
               WRITE (UNIT=output_unit, FMT="(T8,2I5,(T25,4F15.6))") &
                  iset, ishell, &
                  (gto_basis_set%gcc(ipgf, ishell, iset), &
                   ipgf=1, gto_basis_set%npgf(iset))
            END DO
         END DO

         WRITE (UNIT=output_unit, FMT="(A5)") "CPHI"
         WRITE (UNIT=output_unit, FMT="(12F10.5)") gto_basis_set%cphi
         WRITE (UNIT=output_unit, FMT="(A1)") "SPHI"
         WRITE (UNIT=output_unit, FMT="(12F10.5)") gto_basis_set%sphi
         WRITE (UNIT=output_unit, FMT="(A1)") "SCON"
         WRITE (UNIT=output_unit, FMT="(12F10.5)") gto_basis_set%scon

      END IF

   END SUBROUTINE dump_gto_basis_set

! **************************************************************************************************
!> \brief ...
!> \param sto_basis_set ...
! **************************************************************************************************
   SUBROUTINE allocate_sto_basis_set(sto_basis_set)

      TYPE(sto_basis_set_type), POINTER                  :: sto_basis_set

      CALL deallocate_sto_basis_set(sto_basis_set)

      ALLOCATE (sto_basis_set)

   END SUBROUTINE allocate_sto_basis_set

! **************************************************************************************************
!> \brief ...
!> \param sto_basis_set ...
! **************************************************************************************************
   SUBROUTINE deallocate_sto_basis_set(sto_basis_set)

      TYPE(sto_basis_set_type), POINTER                  :: sto_basis_set

! -------------------------------------------------------------------------

      IF (ASSOCIATED(sto_basis_set)) THEN
         IF (ASSOCIATED(sto_basis_set%symbol)) THEN
            DEALLOCATE (sto_basis_set%symbol)
         END IF
         IF (ASSOCIATED(sto_basis_set%nq)) THEN
            DEALLOCATE (sto_basis_set%nq)
         END IF
         IF (ASSOCIATED(sto_basis_set%lq)) THEN
            DEALLOCATE (sto_basis_set%lq)
         END IF
         IF (ASSOCIATED(sto_basis_set%zet)) THEN
            DEALLOCATE (sto_basis_set%zet)
         END IF

         DEALLOCATE (sto_basis_set)
      END IF
   END SUBROUTINE deallocate_sto_basis_set

! **************************************************************************************************
!> \brief ...
!> \param sto_basis_set ...
!> \param name ...
!> \param nshell ...
!> \param symbol ...
!> \param nq ...
!> \param lq ...
!> \param zet ...
!> \param maxlq ...
!> \param numsto ...
! **************************************************************************************************
   SUBROUTINE get_sto_basis_set(sto_basis_set, name, nshell, symbol, nq, lq, zet, maxlq, numsto)

      TYPE(sto_basis_set_type), INTENT(IN)               :: sto_basis_set
      CHARACTER(LEN=default_string_length), &
         INTENT(OUT), OPTIONAL                           :: name
      INTEGER, INTENT(OUT), OPTIONAL                     :: nshell
      CHARACTER(LEN=6), DIMENSION(:), OPTIONAL, POINTER  :: symbol
      INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: nq, lq
      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: zet
      INTEGER, INTENT(OUT), OPTIONAL                     :: maxlq, numsto

      INTEGER                                            :: iset

      IF (PRESENT(name)) name = sto_basis_set%name
      IF (PRESENT(nshell)) nshell = sto_basis_set%nshell
      IF (PRESENT(symbol)) symbol => sto_basis_set%symbol
      IF (PRESENT(nq)) nq => sto_basis_set%nq
      IF (PRESENT(lq)) lq => sto_basis_set%lq
      IF (PRESENT(zet)) zet => sto_basis_set%zet
      IF (PRESENT(maxlq)) THEN
         maxlq = MAXVAL(sto_basis_set%lq(1:sto_basis_set%nshell))
      END IF
      IF (PRESENT(numsto)) THEN
         numsto = 0
         DO iset = 1, sto_basis_set%nshell
            numsto = numsto + 2*sto_basis_set%lq(iset) + 1
         END DO
      END IF

   END SUBROUTINE get_sto_basis_set

! **************************************************************************************************
!> \brief ...
!> \param sto_basis_set ...
!> \param name ...
!> \param nshell ...
!> \param symbol ...
!> \param nq ...
!> \param lq ...
!> \param zet ...
! **************************************************************************************************
   SUBROUTINE set_sto_basis_set(sto_basis_set, name, nshell, symbol, nq, lq, zet)

      TYPE(sto_basis_set_type), INTENT(INOUT)            :: sto_basis_set
      CHARACTER(LEN=default_string_length), INTENT(IN), &
         OPTIONAL                                        :: name
      INTEGER, INTENT(IN), OPTIONAL                      :: nshell
      CHARACTER(LEN=6), DIMENSION(:), OPTIONAL, POINTER  :: symbol
      INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: nq, lq
      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: zet

      INTEGER                                            :: ns

      IF (PRESENT(name)) sto_basis_set%name = name
      IF (PRESENT(nshell)) sto_basis_set%nshell = nshell
      IF (PRESENT(symbol)) THEN
         ns = SIZE(symbol)
         IF (ASSOCIATED(sto_basis_set%symbol)) DEALLOCATE (sto_basis_set%symbol)
         ALLOCATE (sto_basis_set%symbol(1:ns))
         sto_basis_set%symbol(:) = symbol(:)
      END IF
      IF (PRESENT(nq)) THEN
         ns = SIZE(nq)
         CALL reallocate(sto_basis_set%nq, 1, ns)
         sto_basis_set%nq = nq(:)
      END IF
      IF (PRESENT(lq)) THEN
         ns = SIZE(lq)
         CALL reallocate(sto_basis_set%lq, 1, ns)
         sto_basis_set%lq = lq(:)
      END IF
      IF (PRESENT(zet)) THEN
         ns = SIZE(zet)
         CALL reallocate(sto_basis_set%zet, 1, ns)
         sto_basis_set%zet = zet(:)
      END IF

   END SUBROUTINE set_sto_basis_set

! **************************************************************************************************
!> \brief ...
!> \param element_symbol ...
!> \param basis_set_name ...
!> \param sto_basis_set ...
!> \param para_env ...
!> \param dft_section ...
! **************************************************************************************************
   SUBROUTINE read_sto_basis_set(element_symbol, basis_set_name, sto_basis_set, para_env, dft_section)

      ! Read a Slater-type orbital (STO) basis set from the database file.

      CHARACTER(LEN=*), INTENT(IN)                       :: element_symbol, basis_set_name
      TYPE(sto_basis_set_type), INTENT(INOUT)            :: sto_basis_set
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(section_vals_type), POINTER                   :: dft_section

      CHARACTER(LEN=10)                                  :: nlsym
      CHARACTER(LEN=2)                                   :: lsym
      CHARACTER(LEN=240)                                 :: line
      CHARACTER(LEN=242)                                 :: line2
      CHARACTER(len=default_path_length)                 :: basis_set_file_name, tmp
      CHARACTER(LEN=default_path_length), DIMENSION(:), &
         POINTER                                         :: cbasis
      CHARACTER(LEN=LEN(basis_set_name))                 :: bsname
      CHARACTER(LEN=LEN(basis_set_name)+2)               :: bsname2
      CHARACTER(LEN=LEN(element_symbol))                 :: symbol
      CHARACTER(LEN=LEN(element_symbol)+2)               :: symbol2
      INTEGER                                            :: ibasis, irep, iset, nbasis, nq, nset, &
                                                            strlen1, strlen2
      LOGICAL                                            :: basis_found, found, match
      REAL(KIND=dp)                                      :: zet
      TYPE(cp_parser_type)                               :: parser

      line = ""
      line2 = ""
      symbol = ""
      symbol2 = ""
      bsname = ""
      bsname2 = ""

      nbasis = 1

      sto_basis_set%name = basis_set_name
      CALL section_vals_val_get(dft_section, "BASIS_SET_FILE_NAME", &
                                n_rep_val=nbasis)
      ALLOCATE (cbasis(nbasis))
      DO ibasis = 1, nbasis
         CALL section_vals_val_get(dft_section, "BASIS_SET_FILE_NAME", &
                                   i_rep_val=ibasis, c_val=cbasis(ibasis))
         basis_set_file_name = cbasis(ibasis)
         tmp = basis_set_file_name
         CALL uppercase(tmp)
      END DO

      ! Search for the requested basis set in the basis set file
      ! until the basis set is found or the end of file is reached

      basis_found = .FALSE.
      basis_loop: DO ibasis = 1, nbasis
         IF (basis_found) EXIT basis_loop
         basis_set_file_name = cbasis(ibasis)
         CALL parser_create(parser, basis_set_file_name, para_env=para_env)

         bsname = basis_set_name
         symbol = element_symbol
         irep = 0

         tmp = basis_set_name
         CALL uppercase(tmp)

         IF (tmp /= "NONE") THEN
            search_loop: DO

               CALL parser_search_string(parser, TRIM(bsname), .TRUE., found, line)
               IF (found) THEN
                  CALL uppercase(symbol)
                  CALL uppercase(bsname)

                  match = .FALSE.
                  CALL uppercase(line)
                  ! Check both the element symbol and the basis set name
                  line2 = " "//line//" "
                  symbol2 = " "//TRIM(symbol)//" "
                  bsname2 = " "//TRIM(bsname)//" "
                  strlen1 = LEN_TRIM(symbol2) + 1
                  strlen2 = LEN_TRIM(bsname2) + 1

                  IF ((INDEX(line2, symbol2(:strlen1)) > 0) .AND. &
                      (INDEX(line2, bsname2(:strlen2)) > 0)) match = .TRUE.
                  IF (match) THEN
                     ! Read the basis set information
                     CALL parser_get_object(parser, nset, newline=.TRUE.)
                     sto_basis_set%nshell = nset

                     CALL reallocate(sto_basis_set%nq, 1, nset)
                     CALL reallocate(sto_basis_set%lq, 1, nset)
                     CALL reallocate(sto_basis_set%zet, 1, nset)
                     ALLOCATE (sto_basis_set%symbol(nset))

                     DO iset = 1, nset
                        CALL parser_get_object(parser, nq, newline=.TRUE.)
                        CALL parser_get_object(parser, lsym)
                        CALL parser_get_object(parser, zet)
                        sto_basis_set%nq(iset) = nq
                        sto_basis_set%zet(iset) = zet
                        WRITE (nlsym, "(I2,A)") nq, TRIM(lsym)
                        sto_basis_set%symbol(iset) = TRIM(nlsym)
                        SELECT CASE (TRIM(lsym))
                        CASE ("S", "s")
                           sto_basis_set%lq(iset) = 0
                        CASE ("P", "p")
                           sto_basis_set%lq(iset) = 1
                        CASE ("D", "d")
                           sto_basis_set%lq(iset) = 2
                        CASE ("F", "f")
                           sto_basis_set%lq(iset) = 3
                        CASE ("G", "g")
                           sto_basis_set%lq(iset) = 4
                        CASE ("H", "h")
                           sto_basis_set%lq(iset) = 5
                        CASE ("I", "i", "J", "j")
                           sto_basis_set%lq(iset) = 6
                        CASE ("K", "k")
                           sto_basis_set%lq(iset) = 7
                        CASE ("L", "l")
                           sto_basis_set%lq(iset) = 8
                        CASE ("M", "m")
                           sto_basis_set%lq(iset) = 9
                        CASE DEFAULT
                           CALL cp_abort(__LOCATION__, &
                                         "The requested basis set <"//TRIM(bsname)// &
                                         "> for element <"//TRIM(symbol)//"> has an invalid component: ")
                        END SELECT
                     END DO

                     basis_found = .TRUE.
                     EXIT search_loop
                  END IF
               ELSE
                  EXIT search_loop
               END IF
            END DO search_loop
         ELSE
            match = .FALSE.
         END IF

         CALL parser_release(parser)

      END DO basis_loop

      IF (tmp /= "NONE") THEN
         IF (.NOT. basis_found) THEN
            basis_set_file_name = ""
            DO ibasis = 1, nbasis
               basis_set_file_name = TRIM(basis_set_file_name)//"<"//TRIM(cbasis(ibasis))//"> "
            END DO
            CALL cp_abort(__LOCATION__, &
                          "The requested basis set <"//TRIM(bsname)// &
                          "> for element <"//TRIM(symbol)//"> was not "// &
                          "found in the basis set files "// &
                          TRIM(basis_set_file_name))
         END IF
      END IF
      DEALLOCATE (cbasis)

   END SUBROUTINE read_sto_basis_set

! **************************************************************************************************
!> \brief ...
!> \param sto_basis_set ...
!> \param gto_basis_set ...
!> \param ngauss ...
!> \param ortho ...
! **************************************************************************************************
   SUBROUTINE create_gto_from_sto_basis(sto_basis_set, gto_basis_set, ngauss, ortho)

      TYPE(sto_basis_set_type), INTENT(IN)               :: sto_basis_set
      TYPE(gto_basis_set_type), POINTER                  :: gto_basis_set
      INTEGER, INTENT(IN), OPTIONAL                      :: ngauss
      LOGICAL, INTENT(IN), OPTIONAL                      :: ortho

      INTEGER, PARAMETER                                 :: maxng = 6

      CHARACTER(LEN=default_string_length)               :: name, sng
      INTEGER                                            :: ipgf, iset, maxl, ng, nset, nshell
      INTEGER, DIMENSION(:), POINTER                     :: lq, nq
      LOGICAL                                            :: do_ortho
      REAL(KIND=dp), DIMENSION(:), POINTER               :: zet
      REAL(KIND=dp), DIMENSION(maxng)                    :: gcc, zetg

      ng = 6
      IF (PRESENT(ngauss)) ng = ngauss
      IF (ng > maxng) CPABORT("Too many Gaussian primitives requested")
      do_ortho = .FALSE.
      IF (PRESENT(ortho)) do_ortho = ortho

      CALL allocate_gto_basis_set(gto_basis_set)

      CALL get_sto_basis_set(sto_basis_set, name=name, nshell=nshell, nq=nq, &
                             lq=lq, zet=zet)

      maxl = MAXVAL(lq)
      CALL init_orbital_pointers(maxl)

      CALL integer_to_string(ng, sng)
      gto_basis_set%name = TRIM(name)//"_STO-"//TRIM(sng)//"G"

      nset = nshell
      gto_basis_set%nset = nset
      CALL reallocate(gto_basis_set%lmax, 1, nset)
      CALL reallocate(gto_basis_set%lmin, 1, nset)
      CALL reallocate(gto_basis_set%npgf, 1, nset)
      CALL reallocate(gto_basis_set%nshell, 1, nset)
      CALL reallocate(gto_basis_set%n, 1, 1, 1, nset)
      CALL reallocate(gto_basis_set%l, 1, 1, 1, nset)
      CALL reallocate(gto_basis_set%zet, 1, ng, 1, nset)
      CALL reallocate(gto_basis_set%gcc, 1, ng, 1, 1, 1, nset)

      DO iset = 1, nset
         CALL get_sto_ng(zet(iset), ng, nq(iset), lq(iset), zetg, gcc)
         gto_basis_set%lmax(iset) = lq(iset)
         gto_basis_set%lmin(iset) = lq(iset)
         gto_basis_set%npgf(iset) = ng
         gto_basis_set%nshell(iset) = 1
         gto_basis_set%n(1, iset) = lq(iset) + 1
         gto_basis_set%l(1, iset) = lq(iset)
         DO ipgf = 1, ng
            gto_basis_set%gcc(ipgf, 1, iset) = gcc(ipgf)
            gto_basis_set%zet(ipgf, iset) = zetg(ipgf)
         END DO
      END DO

      CALL process_gto_basis(gto_basis_set, do_ortho, nset, maxl)

   END SUBROUTINE create_gto_from_sto_basis

! **************************************************************************************************
!> \brief ...
!> \param gto_basis_set ...
!> \param do_ortho ...
!> \param nset ...
!> \param maxl ...
! **************************************************************************************************
   SUBROUTINE process_gto_basis(gto_basis_set, do_ortho, nset, maxl)

      TYPE(gto_basis_set_type), POINTER                  :: gto_basis_set
      LOGICAL, INTENT(IN), OPTIONAL                      :: do_ortho
      INTEGER, INTENT(IN)                                :: nset, maxl

      INTEGER                                            :: i1, i2, ico, iset, jset, l, lshell, m, &
                                                            maxco, ncgf, ng, ngs, np, nsgf
      INTEGER, DIMENSION(0:10)                           :: mxf
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: gal, zal, zll

      ng = gto_basis_set%npgf(1)
      DO iset = 1, nset
         IF ((ng /= gto_basis_set%npgf(iset)) .AND. do_ortho) &
            CPABORT("different number of primitves")
      END DO

      IF (do_ortho) THEN
         mxf = 0
         DO iset = 1, nset
            l = gto_basis_set%l(1, iset)
            mxf(l) = mxf(l) + 1
         END DO
         m = MAXVAL(mxf)
         IF (m > 1) THEN
            ALLOCATE (gal(ng, nset), zal(ng, nset), zll(m*ng, 0:maxl))
            DO iset = 1, nset
               zal(1:ng, iset) = gto_basis_set%zet(1:ng, iset)
               gal(1:ng, iset) = gto_basis_set%gcc(1:ng, 1, iset)
            END DO
            CALL reallocate(gto_basis_set%zet, 1, m*ng, 1, nset)
            CALL reallocate(gto_basis_set%gcc, 1, m*ng, 1, 1, 1, nset)
            DO iset = 1, nset
               l = gto_basis_set%l(1, iset)
               gto_basis_set%npgf(iset) = ng*mxf(l)
            END DO
            gto_basis_set%zet = 0.0_dp
            gto_basis_set%gcc = 0.0_dp
            zll = 0.0_dp
            mxf = 0
            DO iset = 1, nset
               l = gto_basis_set%l(1, iset)
               mxf(l) = mxf(l) + 1
               i1 = mxf(l)*ng - ng + 1
               i2 = mxf(l)*ng
               zll(i1:i2, l) = zal(1:ng, iset)
               gto_basis_set%gcc(i1:i2, 1, iset) = gal(1:ng, iset)
            END DO
            DO iset = 1, nset
               l = gto_basis_set%l(1, iset)
               gto_basis_set%zet(:, iset) = zll(:, l)
            END DO
            DO iset = 1, nset
               l = gto_basis_set%l(1, iset)
               DO jset = 1, iset - 1
                  IF (gto_basis_set%l(1, iset) == l) THEN
                     m = mxf(l)*ng
                     CALL orthofun(gto_basis_set%zet(1:m, iset), gto_basis_set%gcc(1:m, 1, iset), &
                                   gto_basis_set%gcc(1:m, 1, jset), l)
                  END IF
               END DO
            END DO
            DEALLOCATE (gal, zal, zll)
         END IF
      END IF

      ngs = MAXVAL(gto_basis_set%npgf(1:nset))
      CALL reallocate(gto_basis_set%set_radius, 1, nset)
      CALL reallocate(gto_basis_set%pgf_radius, 1, ngs, 1, nset)
      CALL reallocate(gto_basis_set%first_cgf, 1, 1, 1, nset)
      CALL reallocate(gto_basis_set%first_sgf, 1, 1, 1, nset)
      CALL reallocate(gto_basis_set%last_cgf, 1, 1, 1, nset)
      CALL reallocate(gto_basis_set%last_sgf, 1, 1, 1, nset)
      CALL reallocate(gto_basis_set%ncgf_set, 1, nset)
      CALL reallocate(gto_basis_set%nsgf_set, 1, nset)

      maxco = 0
      ncgf = 0
      nsgf = 0

      DO iset = 1, nset
         gto_basis_set%ncgf_set(iset) = 0
         gto_basis_set%nsgf_set(iset) = 0
         lshell = gto_basis_set%l(1, iset)
         gto_basis_set%first_cgf(1, iset) = ncgf + 1
         ncgf = ncgf + nco(lshell)
         gto_basis_set%last_cgf(1, iset) = ncgf
         gto_basis_set%ncgf_set(iset) = &
            gto_basis_set%ncgf_set(iset) + nco(lshell)
         gto_basis_set%first_sgf(1, iset) = nsgf + 1
         nsgf = nsgf + nso(lshell)
         gto_basis_set%last_sgf(1, iset) = nsgf
         gto_basis_set%nsgf_set(iset) = &
            gto_basis_set%nsgf_set(iset) + nso(lshell)
         ngs = gto_basis_set%npgf(iset)
         maxco = MAX(maxco, ngs*ncoset(lshell))
      END DO

      gto_basis_set%ncgf = ncgf
      gto_basis_set%nsgf = nsgf

      CALL reallocate(gto_basis_set%cphi, 1, maxco, 1, ncgf)
      CALL reallocate(gto_basis_set%sphi, 1, maxco, 1, nsgf)
      CALL reallocate(gto_basis_set%scon, 1, maxco, 1, nsgf)
      CALL reallocate(gto_basis_set%lx, 1, ncgf)
      CALL reallocate(gto_basis_set%ly, 1, ncgf)
      CALL reallocate(gto_basis_set%lz, 1, ncgf)
      CALL reallocate(gto_basis_set%m, 1, nsgf)
      CALL reallocate(gto_basis_set%norm_cgf, 1, ncgf)
      ALLOCATE (gto_basis_set%cgf_symbol(ncgf))
      ALLOCATE (gto_basis_set%sgf_symbol(nsgf))

      ncgf = 0
      nsgf = 0

      DO iset = 1, nset
         lshell = gto_basis_set%l(1, iset)
         np = lshell + 1
         DO ico = ncoset(lshell - 1) + 1, ncoset(lshell)
            ncgf = ncgf + 1
            gto_basis_set%lx(ncgf) = indco(1, ico)
            gto_basis_set%ly(ncgf) = indco(2, ico)
            gto_basis_set%lz(ncgf) = indco(3, ico)
            gto_basis_set%cgf_symbol(ncgf) = &
               cgf_symbol(np, [gto_basis_set%lx(ncgf), &
                               gto_basis_set%ly(ncgf), &
                               gto_basis_set%lz(ncgf)])
         END DO
         DO m = -lshell, lshell
            nsgf = nsgf + 1
            gto_basis_set%m(nsgf) = m
            gto_basis_set%sgf_symbol(nsgf) = sgf_symbol(np, lshell, m)
         END DO
      END DO

      gto_basis_set%norm_type = -1

   END SUBROUTINE process_gto_basis

! **************************************************************************************************
!> \brief ...
!> \param zet ...
!> \param co ...
!> \param cr ...
!> \param l ...
! **************************************************************************************************
   SUBROUTINE orthofun(zet, co, cr, l)
      REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: zet
      REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: co, cr
      INTEGER, INTENT(IN)                                :: l

      REAL(KIND=dp)                                      :: ss

      CALL aovlp(l, zet, cr, cr, ss)
      cr(:) = cr(:)/SQRT(ss)
      CALL aovlp(l, zet, co, cr, ss)
      co(:) = co(:) - ss*cr(:)
      CALL aovlp(l, zet, co, co, ss)
      co(:) = co(:)/SQRT(ss)

   END SUBROUTINE orthofun

! **************************************************************************************************
!> \brief ...
!> \param l ...
!> \param zet ...
!> \param ca ...
!> \param cb ...
!> \param ss ...
! **************************************************************************************************
   SUBROUTINE aovlp(l, zet, ca, cb, ss)
      INTEGER, INTENT(IN)                                :: l
      REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: zet, ca, cb
      REAL(KIND=dp), INTENT(OUT)                         :: ss

      INTEGER                                            :: i, j, m
      REAL(KIND=dp)                                      :: ab, ai, aj, s00, sss

!
! use init_norm_cgf_orb
!
      m = SIZE(zet)
      ss = 0.0_dp
      DO i = 1, m
         ai = (2.0_dp*zet(i)/pi)**0.75_dp
         DO j = 1, m
            aj = (2.0_dp*zet(j)/pi)**0.75_dp
            ab = 1._dp/(zet(i) + zet(j))
            s00 = ai*aj*(pi*ab)**1.50_dp
            IF (l == 0) THEN
               sss = s00
            ELSEIF (l == 1) THEN
               sss = s00*ab*0.5_dp
            ELSE
               CPABORT("aovlp lvalue")
            END IF
            ss = ss + sss*ca(i)*cb(j)
         END DO
      END DO

   END SUBROUTINE aovlp

! **************************************************************************************************
!> \brief ...
!> \param z ...
!> \param ne ...
!> \param n ...
!> \param l ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION srules(z, ne, n, l)
      ! Slater rules
      INTEGER, INTENT(IN)                                :: z
      INTEGER, DIMENSION(:, :), INTENT(IN)               :: ne
      INTEGER, INTENT(IN)                                :: n, l
      REAL(dp)                                           :: srules

      REAL(dp), DIMENSION(7), PARAMETER :: &
         xns = [1.0_dp, 2.0_dp, 3.0_dp, 3.7_dp, 4.0_dp, 4.2_dp, 4.4_dp]

      INTEGER                                            :: i, l1, l2, m, m1, m2, nn
      REAL(dp)                                           :: s

      s = 0.0_dp
      ! The complete shell
      l1 = MIN(l + 1, 4)
      nn = MIN(n, 7)
      IF (l1 == 1) l2 = 2
      IF (l1 == 2) l2 = 1
      IF (l1 == 3) l2 = 4
      IF (l1 == 4) l2 = 3
      ! Rule a) no contribution from shells further out
      ! Rule b) 0.35 (1s 0.3) from each other electron in the same shell
      IF (n == 1) THEN
         m = ne(1, 1)
         s = s + 0.3_dp*REAL(m - 1, dp)
      ELSE
         m = ne(l1, nn) + ne(l2, nn)
         s = s + 0.35_dp*REAL(m - 1, dp)
      END IF
      ! Rule c) if (s,p) shell 0.85 from each electron with n-1, and 1.0
      !      from all electrons further in
      IF (l1 + l2 == 3) THEN
         IF (nn > 1) THEN
            m1 = ne(1, nn - 1) + ne(2, nn - 1) + ne(3, nn - 1) + ne(4, nn - 1)
            m2 = 0
            DO i = 1, nn - 2
               m2 = m2 + ne(1, i) + ne(2, i) + ne(3, i) + ne(4, I)
            END DO
            s = s + 0.85_dp*REAL(m1, dp) + 1._dp*REAL(m2, dp)
         END IF
      ELSE
         ! Rule d) if (d,f) shell 1.0 from each electron inside
         m = 0
         DO i = 1, nn - 1
            m = m + ne(1, i) + ne(2, i) + ne(3, i) + ne(4, i)
         END DO
         s = s + 1._dp*REAL(m, dp)
      END IF
      ! Slater exponent is (Z-S)/NS
      srules = (REAL(z, dp) - s)/xns(nn)
   END FUNCTION srules

! **************************************************************************************************
!> \brief sort basis sets w.r.t. radius
!> \param basis_set ...
!> \param sort_method ...
! **************************************************************************************************
   SUBROUTINE sort_gto_basis_set(basis_set, sort_method)
      TYPE(gto_basis_set_type), INTENT(INOUT)            :: basis_set
      INTEGER, INTENT(IN)                                :: sort_method

      CHARACTER(LEN=12), DIMENSION(:), POINTER           :: cgf_symbol
      CHARACTER(LEN=6), DIMENSION(:), POINTER            :: sgf_symbol
      INTEGER :: ic, ic_max, icgf, icgf_new, icgf_old, ico, is, is_max, iset, isgf, isgf_new, &
         isgf_old, ishell, lshell, maxco, maxpgf, maxshell, mm, nc, ncgf, ns, nset
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: sort_index
      INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: icgf_set, isgf_set
      INTEGER, DIMENSION(:), POINTER                     :: lx, ly, lz, m, npgf
      REAL(dp), ALLOCATABLE, DIMENSION(:)                :: tmp
      REAL(dp), DIMENSION(:), POINTER                    :: set_radius
      REAL(dp), DIMENSION(:, :), POINTER                 :: zet
      REAL(KIND=dp), DIMENSION(:), POINTER               :: norm_cgf
      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: cphi, scon, sphi

      NULLIFY (set_radius, zet)

      IF (sort_method == basis_sort_default) RETURN

      CALL get_gto_basis_set(gto_basis_set=basis_set, &
                             nset=nset, &
                             maxshell=maxshell, &
                             maxpgf=maxpgf, &
                             maxco=maxco, &
                             ncgf=ncgf, &
                             npgf=npgf, &
                             set_radius=set_radius, &
                             zet=zet)

      ALLOCATE (sort_index(nset))
      ALLOCATE (tmp(nset))
      SELECT CASE (sort_method)
      CASE (basis_sort_zet)
         DO iset = 1, nset
            tmp(iset) = MINVAL(basis_set%zet(:npgf(iset), iset))
         END DO
      CASE DEFAULT
         CPABORT("Request basis sort criterion not implemented.")
      END SELECT

      CALL sort(tmp(1:nset), nset, sort_index)

      ic_max = 0
      is_max = 0
      DO iset = 1, nset
         ic = 0
         is = 0
         DO ishell = 1, basis_set%nshell(iset)
            DO ico = 1, nco(basis_set%l(ishell, iset))
               ic = ic + 1
               IF (ic > ic_max) ic_max = ic
            END DO
            lshell = basis_set%l(ishell, iset)
            DO mm = -lshell, lshell
               is = is + 1
               IF (is > is_max) is_max = is
            END DO
         END DO
      END DO

      icgf = 0
      isgf = 0
      ALLOCATE (icgf_set(nset, ic_max))
      icgf_set(:, :) = 0
      ALLOCATE (isgf_set(nset, is_max))
      isgf_set(:, :) = 0

      DO iset = 1, nset
         ic = 0
         is = 0
         DO ishell = 1, basis_set%nshell(iset)
            DO ico = 1, nco(basis_set%l(ishell, iset))
               icgf = icgf + 1
               ic = ic + 1
               icgf_set(iset, ic) = icgf
            END DO
            lshell = basis_set%l(ishell, iset)
            DO mm = -lshell, lshell
               isgf = isgf + 1
               is = is + 1
               isgf_set(iset, is) = isgf
            END DO
         END DO
      END DO

      ALLOCATE (cgf_symbol(SIZE(basis_set%cgf_symbol)))
      ALLOCATE (norm_cgf(SIZE(basis_set%norm_cgf)))
      ALLOCATE (lx(SIZE(basis_set%lx)))
      ALLOCATE (ly(SIZE(basis_set%ly)))
      ALLOCATE (lz(SIZE(basis_set%lz)))
      ALLOCATE (cphi(SIZE(basis_set%cphi, 1), SIZE(basis_set%cphi, 2)))
      cphi = 0.0_dp
      ALLOCATE (sphi(SIZE(basis_set%sphi, 1), SIZE(basis_set%sphi, 2)))
      sphi = 0.0_dp
      ALLOCATE (scon(SIZE(basis_set%scon, 1), SIZE(basis_set%scon, 2)))
      scon = 0.0_dp

      ALLOCATE (sgf_symbol(SIZE(basis_set%sgf_symbol)))
      ALLOCATE (m(SIZE(basis_set%m)))

      icgf_new = 0
      isgf_new = 0
      DO iset = 1, nset
         DO ic = 1, ic_max
            icgf_old = icgf_set(sort_index(iset), ic)
            IF (icgf_old == 0) CYCLE
            icgf_new = icgf_new + 1
            norm_cgf(icgf_new) = basis_set%norm_cgf(icgf_old)
            lx(icgf_new) = basis_set%lx(icgf_old)
            ly(icgf_new) = basis_set%ly(icgf_old)
            lz(icgf_new) = basis_set%lz(icgf_old)
            cphi(:, icgf_new) = basis_set%cphi(:, icgf_old)
            cgf_symbol(icgf_new) = basis_set%cgf_symbol(icgf_old)
         END DO
         DO is = 1, is_max
            isgf_old = isgf_set(sort_index(iset), is)
            IF (isgf_old == 0) CYCLE
            isgf_new = isgf_new + 1
            m(isgf_new) = basis_set%m(isgf_old)
            sphi(:, isgf_new) = basis_set%sphi(:, isgf_old)
            scon(:, isgf_new) = basis_set%scon(:, isgf_old)
            sgf_symbol(isgf_new) = basis_set%sgf_symbol(isgf_old)
         END DO
      END DO

      DEALLOCATE (basis_set%cgf_symbol)
      basis_set%cgf_symbol => cgf_symbol
      DEALLOCATE (basis_set%norm_cgf)
      basis_set%norm_cgf => norm_cgf
      DEALLOCATE (basis_set%lx)
      basis_set%lx => lx
      DEALLOCATE (basis_set%ly)
      basis_set%ly => ly
      DEALLOCATE (basis_set%lz)
      basis_set%lz => lz
      DEALLOCATE (basis_set%cphi)
      basis_set%cphi => cphi
      DEALLOCATE (basis_set%sphi)
      basis_set%sphi => sphi
      DEALLOCATE (basis_set%scon)
      basis_set%scon => scon

      DEALLOCATE (basis_set%m)
      basis_set%m => m
      DEALLOCATE (basis_set%sgf_symbol)
      basis_set%sgf_symbol => sgf_symbol

      basis_set%lmax = basis_set%lmax(sort_index)
      basis_set%lmin = basis_set%lmin(sort_index)
      basis_set%npgf = basis_set%npgf(sort_index)
      basis_set%nshell = basis_set%nshell(sort_index)
      basis_set%ncgf_set = basis_set%ncgf_set(sort_index)
      basis_set%nsgf_set = basis_set%nsgf_set(sort_index)

      basis_set%n(:, :) = basis_set%n(:, sort_index)
      basis_set%l(:, :) = basis_set%l(:, sort_index)
      basis_set%zet(:, :) = basis_set%zet(:, sort_index)

      basis_set%gcc(:, :, :) = basis_set%gcc(:, :, sort_index)
      basis_set%set_radius(:) = basis_set%set_radius(sort_index)
      basis_set%pgf_radius(:, :) = basis_set%pgf_radius(:, sort_index)

      nc = 0
      ns = 0
      DO iset = 1, nset
         DO ishell = 1, basis_set%nshell(iset)
            lshell = basis_set%l(ishell, iset)
            basis_set%first_cgf(ishell, iset) = nc + 1
            nc = nc + nco(lshell)
            basis_set%last_cgf(ishell, iset) = nc
            basis_set%first_sgf(ishell, iset) = ns + 1
            ns = ns + nso(lshell)
            basis_set%last_sgf(ishell, iset) = ns
         END DO
      END DO

   END SUBROUTINE sort_gto_basis_set

END MODULE basis_set_types
